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ABSTRACT 


A  three-dimensional  numerical  model  of  the  Jet  Engine  Hush  House  located  at  Na¬ 
val  Air  Station  Jacksonville.  Florida  was  developed  using  the  PHOENICS  computer 
code.  The  PHOENICS  code  was  used  to  determine  the  steady  state  aerothermal  char¬ 
acteristics  in  the  Hush  House  produced  by  a  Na\'y  F-4  (Phantom  11)  J-79-GE-S  gas 
turbine  engine  with  afterburner.  The  PHOENICS  generated  pressure,  velocity,  and 
temperature  fields  of  the  test  cell  facility  based  on  the  k-£  turbulence  model  are  com¬ 
pared  with  results  generated  by  a  model  incorporating  a  temperature  dependent 
kinematic  turbulent  viscosity.  In  addition,  for  the  k-c  turbulence  model,  the  effect  of 
augmenter  tube  length  on  sytem  behavior  was  investigated.  Lessons  learned  and  rec- 
onunendatiens  in  obtaining  a  converged  solution  are  included  to  aid  in  further  applica¬ 
tions  of  the  PHOENICS  code. 
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I.  INTRODUCTION 


Todays  sophisticated  aircraft  engines,  both  new  and  overhauled,  require  f”'!  opera¬ 
tional  testing  prior  to  being  placed  in  ser  .ice.  Jet  engine  test  facilities  are  utilized  to 
safely  and  eOectisely  monitor  aircraft  and  iet  engine  performance  characteristics.  In 
addition  to  providing  an  environment  for  acquiring  test  data,  the  te.>t  cells  supply  cli  an. 
distortion-free  airflow  to  the  engine  inlet,  decrease  noise  pollution  by  decreasing  the 
kinetic  energy  of  the  jet,  and  cool  the  exhaust  prior  to  exiting  to  the  atmosphere  [Rel'. 

1;  p,  1], 

There  are  five  types  of  tests  facilities  currently  used  by  the  N.AV^’. 

1.  The  Run-up  Pad--wherc  the  aircraft  is  tied  down  in  an  unprotected  environment 
to  conauct  power  checks. 

2.  Sound  Suppressors--equipment  packages  that  reduce  the  noise  of  aircraft  on  the 
Run-Up  Pad. 

3.  Open  Test  Stand  (OTSb-outdoor  engine  test  stand  and  portable  sensing  equip¬ 
ment. 

4.  Hush  Plouse-a  total  aircraft  enclosure  capable  of  around  the  clock  testing. 

5.  Fully  Enclosed  Test  Cell(TC)”an  enclosed  facility  for  complex  testing  of  engines 
out  of  the  airframe. 

The  run-up  pad.  sound  suppressc..  and  the  open  test  stand  facilities  although  relatively 
inexpensive  are  highly  inefllcient  and  do  not  support  military  readiness  needs.  Noise  and 
air  pollution  factors  determine  the  allowable  locations  of  the  facilities  and  the  environ¬ 
mental  conditions  control  the  time  of  operation  of  some  of  the  facilities.  The  hush  house 
and  the  fully  enclosed  test  cell  can  operate  independent  from  environmental  conditions, 
and  since  noise  is  reduced,  around  the  clock  testing  is  allowed.  [Ref  2:  p.  2] 

As  of  .March  1984.  75  Navy  and  .Marine  activities  have  76  fully  enclosed  test  cells 
and  fis  e  hush  houses.  Only  60  out  of  the  76  test  cells  are  fully  operational.  Most  of  the 
test  cells  are  over  27  years  old.  .Modern  jet  engines  operate  at  much  higher  temperatures. 
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mass  flows,  and  velocities  than  they  did  20  years  ago  when  most  of  the  test  cells  were 
designed  and  built.  Due  to  the  age  of  the  facilities  and  high  gas  temperatures  and  water 
impingement,  cracking  and  spalling  of  the  test  cell  concrete  walls  continues  to  occur. 
In  addition  to  the  structural  degradation  of  the  test  facilities  many  of  the  test  cells  are 
obsolete  in  that  the\  are  unable  to  support  the  increased  size  and  wide  variety  of  the 
engines  requiring  testing.  [Ref  2:  pp.  3-4] 

Noise,  air.  and  water  pollution  control  at  test  facilities  is  under  increasing  scrutiny. 
Prior  to  1976,  military  test  cells  were  exempt  from  state  and  local  regulations  pursuant 
to  section  233  of  the  Clean  Air  Act.  The  Federal  Court  decided  in  1981,  that  jest  engine 
test  cells  ucre  to  adhere  to  the  requirements  of  the  Clean  Air  Act.  The  absence  and 
degradation  of  exsisting  pollution  control  systems  has  increased  the  risk  of  physical  in¬ 
jury  to  the  test  cell  facility  personnel  as  well  as  damage  risk  to  test  cells,  equipment,  and 
aircraft.  In  1981.  the  Navy  Lnvironmental  Health  Center  cited  hearing  loss  and  respir¬ 
atory  damage  for  Slid  million  with  an  anticipated  claim  growth  of  9'^o  each  year. 
Fquipment  and  aircraft  arc  subjected  to  acoustic  fatigue  which  could  possibly  shorten 
operating  life.  [Ref  2:  p.  3] 

It  is  the  Navy's  goal  to  provide  standardized  facilities  for  testing  aircraft  and  aircraft 
jet  engines  that  will  incorporate  the  latest  technology  and  which  will  adhere  to  all  local, 
state,  and  federal  regulations  concerning  environmental  safety.  The  current  jet  engine 
rest  cell  facilities  were  designed  empirically,  based  upon  the  designer's  experience  and 
occasional  aerothcrmal  data  obtained  from  existing  test  cells.  The  test  cell  facilities  in 
general  performed  well  for  the  type  engines  that  they  were  designed  for.  Due  to  the  in¬ 
crease  in  conrplexity  and  power  of  modern  jet  engines  and  aircraft  the  test  cell  facility 
design  procedures  of  the  past  are  inadequate.  The  design  of  a  modern  jet  engine  test 
facility  must  incorporate  the  eflccts  of  multiple  engines  producing  nonconcentric  flows, 
turbulence,  supersonic  velocities,  extreme  temperature  gradients,  and  awkward 


geometries.  In  addition,  the  new  test  cells  must  be  able  to  accomodate  a  wide  varictv 
of  engines  and  aircraft  each  with  its  own  characteristic  exhaust  temperature,  exit 
velocity,  and  mass  flow  rate  [Ref  1;  p.  1). 

The  design  of  modern  test  cell  facilities  requires  a  more  sophisticated  design  proce¬ 
dure.  A  numerical  simulation  of  the  test  cell  can  provide  the  design  engineer  with  critical 
performance  data  from  existing  test  cells  to  aid  in  the  design  of  new  test  cell  facilities. 
The  numerical  model  provides  a  cost  and  time  eflicient  means  of  obtaining  the  per¬ 
formance  data  of  existing  test  cell  facility  designs,  and  this  information  can  be  incorpo¬ 
rated  in  the  design  of  new  test  facilities. 

The  PHOENICS  Code  is  used  to  numerically  model  the  aerotherrnal  characteristics 
in  the  Hush  House  Facility,  located  at  Naval  Air  Station  Jacksonville.  Florida,  due  to 
the  operation  of  a  Nacy  j-79-GF-8  gas  turbine  engine. 


II.  PROBLEM  DESCRIPTION 

A.  HUSH  HOUSE  CHARACTERIZATION 

In  1980  the  Air  Force  purchased  79  standardized  demountable  Hush  Houses  with 
test  cell  capabilities.  The  large  scale  purchase,  standardized  design,  and  simple 
contruction  resulted  in  a  relatively  low  cost  of  S2.2  million  per  Hush  House.  The  Hush 
House  located  at  Naval  Air  Station  Jacksonville,  Florida  is  of  the  same  design  as  the 
Air  Force  Hush  House  except  for  the  additional  AFFF  Firefighting  system  installed  by 
the  Navy.  Neither  the  .Air  Force's  nor  the  Navy's  Hush  House  design  incorporates  air 
pollution  control  systems.  Although  the  cost  of  the  initial  large  scale  purchase  by  the 
.Air  Force  was  low,  the  overall  long  term  cost  was  significantly  higher  due  to  high 
maintenance  repair  costs.  It  appears  that  the  purchase  was  made  in  haste  and  not  fully 
researched.  The  Nasw  could  benefit  from  a  similar  large  scale  purchase  providing  that 
a  high  performance  low  cost  maintenance  system  can  be  designed.  [Ref  2:  P.  1] 

The  basic  operation  of  the  Hush  House  consists  of  four  major  sections:  intakes,  test 
bay,  augmenter  tube,  and  exhaust  stack  as  shown  in  Figure  1  on  page  6  and  Figure  2 
on  page  7.  The  primarv’  intakes  provide  air  for  the  gas  turbine  engine  intake.  1  he  sec¬ 
ondary  intakes  provide  air  to  mix  with  the  engine  exhaust.  The  mixing  of  the  low  tem¬ 
perature  low  velocity  secondaiy  air  with  the  high  temperature  high  velocity  exhaust 
results  in  cooling  the  hot  exhaust  gas  and  reducing  the  exhaust  gas  velocity.  The  test 
bay  is  a  large  unobstructed  area  where  the  entire  aircraft  or  an  engine  mounted  in  a  test 
stand  is  positioned  for  testing.  The  large  test  bay  allows  intake  air  to  smooth  out  into 
a  uniform  profile  prior  to  the  engine  intake  thereby  approximating  true  operating  con¬ 
ditions.  The  augmenter  tube  is  lined  with  acoustic  absorption  material  to  aid  in  acoustic 
energy  dissipation  from  the  engine  and  jet  flow.  The  augmenter  tube  aids  in  the  mixing 
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of  the  engine  exhaust  and  secondan.'  intake  cooling  air.  The  exhaust  stack  changes  the 
flow  direction  from  horizontal  to  vertical  as  the  exhaust  gas  flow  exits  to  the  atmos¬ 
phere.  [Ref  1;  p.  2] 


The  overall  dimensions  of  the  test  bay  are  2S.7  meters  (94  feet  )  wide  by  25.5  meters 
(83.7  feet)  long  by  8.6  meters  (28.2  feet)  high.  The  24.4  meter  (76  feet)  long  oval  shaped 
augmenter  tube  is  centered  at  the  rear  end  of  the  test  bay  with  the  floor  of  the  tube  at 
the  same  level  as  the  test  bay  floor.  The  oval  augmenter  tube  cross  section  is  roughly 
4.9  meters  (16  feet)  wide  by  4.0  meters  (13.2  feet)  high.  The  dimensions  are  constant 
throughout  the  tube  length.  The  exhaust  stack  at  the  end  of  the  augmenter  tube  is 
simply  a  45  degree  ramp  ending  just  short  of  the  top  of  the  augmenter  tube.  The  exhaust 
opening  is  a  8.2  meter  (27  feet)  by  8.2  meter  (27  feet)  horizontal  plane.  The  aircraft  is 
positioned  in  the  test  bay  with  the  engine  outlet  approximately  4.6  meters  (15  feet)  from 
the  augmenter  tube  entrance  (see  Figure  1  on  page  6  and  Figure  2  on  page  7).  General 
operating  procedures  for  conducting  engine  tests  are  described  by  Nicolaus  (1988)  [Ref 
3:  pp.  5-6). 

B.  ISSUES  CONCERNING  THE  HUSH  HOUSE 

The  Air  Force  type  Hush  House  in  general  is  performing  satisfactorly.  Some  struc¬ 
tural  and  operational  problems  do  exist  in  the  Hush  House.  The  stainless  steel  corru¬ 
gated  sheets  that  overlay  the  acoustic  absorption  material  is  subject  to  cracking  along 
the  welds.  During  long  periods  of  idle  engine  operation  a  portion  of  the  exhaust  gas  flow 
is  deflected  away  from  the  augmenter  tube  entrance.  The  exhaust  gas  rises  to  the  top 
of  the  test  bay  activating  the  high  heat  sensors  [Ref  3:  p.  7].  Vibration  of  the  entire 
Hush  House  structure  is  present  during  engine  testing.  It  is  not  exactly  clear  what  effect 
the  vibration  has  on  the  structural  integrity  of  the  Hush  Houses.  Some  of  the  Hush 
Houses  require  inspection  and  retightening  of  the  bolts  that  fasten  the  acoustic  absorp¬ 
tion  material  to  the  walls  of  the  augmenter  tube  [Ref  2:  p.  4). 
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III.  THE  PHOENICS  CODE 


A.  GENERAL  DESCRIPTION 

PHOENICS  is  a  commercially  available  general  purpose  computer  program  copy 
righted  by  CHAM  (Concentration  Heat  and  Momentum)  Limited  which  was  founded 
in  1972  by  Prof.  D.  Brian  Spalding  of  Imperial  College  in  London,  England. 
PHOENICS,  an  acronym  for  Parabolic  Hyperbolic  Or  Elliptic  Numerical  Integration 
Code  Series,  simulates  fluid  flotv,  heat  transfer,  and  chetnical  reactions.fRef.  p.  l.Ij 
The  mathematical  basis  of  PHOENICS  is  the  iterative  solution  of  descretized  con¬ 
servation  equations  of  mass,  chemical  species,  momentum,  energy,  turbulence  quantities, 
and  other  related  lluid  properties.  The  general  form  of  the  conservation  equation  (before 
descrctization)  for  which  PHOENICS  solves  can  be  expressed  in  the  form  [Ref  4;  p,  6.2): 

-Tj  (p<t>)  +  dh(/}  J'(p-r =  S^  (3.1) 


where: 

(/)  =  dependent  variable 

=  exchange  coefficient  for  </> 
s  source  of  4>  per  unit  volume 
P  =  density 

I '  =  velocity  vector 

For  example,  in  the  energy  conservation  equation  the  dependent  variable  4)  would  be 
enthalpy  and  the  exchange  coefficient  E*.  would  be  the  effective  thermal  conductivty  di¬ 
vided  by  the  constant  pressure  specific  heat.  The  partial  differential  equations(p.d.e.) 
are  coupled  and  nonlinear.  The  coupled,  nonlinear  conservation  p.d.e's  can  only  be 
solved  (except  for  some  very  simple  cases  )  by  approximate  numerical  methods.  Typical 
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methods  for  their  approximate  numerical  solution  includes  the  finite  difference  and  finite 
element  methods.  In  the  present  case,  the  control  volume  method  has  features  of  both 
the  finite  difference  and  finite  element  methods.  All  these  methods  transform  a  partial 
differential  equation  into  system  of  algebraic  equations  as  follows: 


Op<j)p  —  a\4>\+  +  U£</)£  +  aii4>l{  +  ^Lf^L  "b  ^7^r"b  ^  ('’•“) 


where: 

a  =  influence  coefficients 

b  =  a  representation  of  the  source  of  <f)  for  the  cell 

The  subscripts  N,  S,  E,  W,  H,  and  L  denote  the  North,  South,  East, West.  High,  and 
Low  neighboring  cell  faces  respectively  as  shown  in  Figure  3  on  page  13.  The  subscript 
P  denotes  the  grid  node  at  a  particular  location  and  the  T  subscript  denotes  the  same 
grid  node  at  an  earlier  time.  The  system  domain  is  partitioned  into  a  grid  network  of 
cells.  This  partition  discretizes  the  continous  domain  into  A'c  cells.  Then  each  p.d.e.  is 
transformed  into  an  algebraic  equation  for  each  cell  as  follows: 

ci\<t)y  +  +  an<f>jy+  +  S  (3 ’’'1 

Uy  +  <^5  +  +  «»•+%  + ^ 

Thus  for  example,  the  energy  p.d.e.  is  transformed  into  A'c  algebraic  equations  which  are 
then  solved  for  the  value  of  enthalpy  (  or  temperature  )  at  the  center  of  each  of  the  .NV 
cells.  It  should  be  noted  that  these  algebraic  equations  are  also  coupled,  and  therefore 
cannot  be  solved  without  solving  the  discretized  algebraic  equations  for  all  of  the  de¬ 
pendent  variables,  The  algebraic  equations  are  numerous  and  dependent  on  similar 
equations  from  neighboring  cells.  The  number  of  algebraic  equations  is  equal  to  the 
number  of  dependent  variables  times  the  number  of  cells,  A^  The  Hush  House  utilizes 
4800  cells.  Solving  for  seven  dependent  variables  in  each  cell  equates  to  the  solution  of 
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33600  equations.  PHOENICS  uses  the  SIMPLE  algorithm  requiring  a  guess  and  cor¬ 
rection  scheme  of  which  the  objective  is  to  reduce  the  imbalance  between  the  left  and 
right  hand  sides  of  the  discretized  equation  to  a  negligible  magnitude  [Ref  5:  p.  2..'']. 
The  specific  methods  of  reducing  the  imbalances  (referred  to  as  residuals)  and  thereby 
obtaining  a  converged  solution  is  discussed  in  Chapter  four. 

The  PHOENICS  computer  code  consists  of  two  main  programs  EAR  TH  and  SAT¬ 
ELLITE.  Figure  4  on  page  14  shows  the  relationships  between  E.\RTH,  SATELLITE, 
and  their  respective  subprograms  and  hardware.  The  EARTH  program  is  the  largest  of 
the  two  and  it  contains  the  coding  for  solving  the  conservation  and  related  physical  law 
equations.  SATELLITE-EARTH  communication  is  uni-directional,  information  can  go 
only  from  SATELLITE  to  EARTH.  Therefore  a  property  which  itself  is  a  function  of 
the  solution,  which  is  achieved  in  EARTH,  cannot  be  accommodated  by  coding  only  in 
SATELLITE.  For  accommodating  solution  dependent  properties,  one  must  use 
GROUND.  The  subroutine  GREXl  is  an  e.xemplary  GROUND  code  that  provides 
common  options  for  general  material  properties  such  as  density  and  viscosity.  Desired 
options  which  are  not  in  GREXl  must  be  programmed  in  GROUND. 

The  SATELLITE  program  transposes  input  data  e.xpressed  in  PHOENICS  input 
language  (PIL)  into  commands  that  will  actvate  EARTH.  A  Q1  file  (quick-input  file, 
number  1 ),  where  the  flow  problem  is  defined,  is  a  formatted  input  file  divided  into  24 
groups; 

GROUP  1  =  Run  title  and  other  preliminaries 
GROUP  2  =  Transcience:  time-step  specification 
GROUP  3  =  X-direction  grid  specification 
GROUP  4  =  Y-direction  grid  specification 
GROUP  5  =  Z-direction  grid  specification 
GROUP  6  =  Body-fitted  coordinates  or  grid  distortion 
GROUP  7  =  Variables  stored,  solved  and  named 
GROUP  8  =  Terms  (in  differential  equations)  and  devices 
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GROUP  9  =  Properties  of  the  medium  {or  media) 

GROUP  10  =  Iterphase-transfer  presses  and  properties 

GROUP  11  s  Initialization  of  variable  or  porosity  fields 

GROUP  12  =  Convection  and  difiusion  adjustments 

GROUP  13  =  Boundary  conditions  and  special  sources 

GROUP  14  =  downstream  pressure  for  PARAB  =  T 

GROUP  15  =  Termination  of  sweeps 

GROUP  16  =  Termination  of  iterations 

GROUP  17  s  Underrelaxation  devices 

GROUP  18  s  Limits  on  variables  or  increments  to  them 

GROUP  19  s  Data  communication  by  SATELLITE  to  GROUND 

GROUP  20  =  Preliminary  printout 

GROUP  21  s  Print  out  of  variables 

GROUP  22  =  Spot-value  printout 

GROUP  23  =  Field  printout  and  plot  control 

GROUP  24  =  Dumps  for  restarts 

A  complete  explanation  of  each  group's  data  structure  and  a  PHOENTCS  Input  Lan¬ 
guage  glossary  can  be  found  in  the  PHOENICS  Beginners  Guide  and  Users  Manual  (sec 
Reference  2). 

B.  MODELING  THE  PROBLEM 

The  main  scope  of  the  thesis  is  a  parametric  study  of  the  Naval  Air  Station, 
Jacksonville,  Florida,  Hush  House  originally  modeled  by  Nicolaus  in  19S8  [Ref.  3].  The 
original  model  of  the  Hush  House  used  the  cartesian  coordinate  grid  system  of  X,  Y,  and 
Z  coordinates.  The  original  three  dimensional  rectangular  grid  consisted  of  40  cells  in 
the  Z  direction  (NZ  =  40,  low  to  high).  12  cells  in  the  Y  direction  (NY=  12,  south  to 
north),  and  10  cells  in  the  X  direction  (NX=  10,  west  to  east).  The  Hush  House  model 
takes  advantage  of  the  Y-Z  plane  of  symmetry’  of  the  actual  Hush  House  and  therefore 
only  half  of  the  Hush  House  is  modeled.  Dividing  the  structure  in  half  reduces  the 
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number  of  equations  to  be  solved  by  half  thereby  reducing  computer  CPU  time.  1  he  te'^t 
bay  is  defined  by  cells  1-19  in  the  Z  direction,  cells  1-12  in  the  direction,  and  cells  1-9 
in  the  X  direction.  The  jet  engine  is  defined  by  cells  12-14  in  the  Z  direction,  cells  4-5 
in  the  Y  direction,  and  cell  1  in  the  X  direction.  The  augmenter  tube  is  defined  by  cells 
lS-27  in  the  Z  direction,  cells  1-8  in  the  Y  direction,  and  cells  1-5  in  the  X  direction.  I  he 
exhaust  stack  is  defined  by  cells  28-39  in  the  Z  direction  ,  cells  1-8  in  the  Y  direction,  and 
cells  1-5  in  the  X  direction.  The  45  degree  ramp  at  the  end  of  the  exhaust  is  modeled 
with  a  step  system.  Figure  5  on  page  15  and  Figure  6  on  page  16  graphically  display 
the  X-Z  and  the  Y-Z  planes  of  the  structure.  (Ref  3:  pp.  25-26) 

The  graphical  results  of  the  temperature  and  velocity  profiles  at  volume  slab  IX  =  1 
for  IY=4  from  the  original  model  shown  in  Figure  7  on  page  17  and  Figure  8  on  page 
IS  indicate  that  the  exhaust  gas  temperature  and  velocity  were  reduced  by  roughly  50 
percent  in  the  first  20  percent  of  the  augmenter  tube  (Ref  3:  pp.  35-36).  The  question 
then  arose  as  to  what  would  be  the  effect  of  a  shorter  augmenter  tube.  The  first 
parametric  study  investigates  the  elTects  of  augmenter  tube  length  on  the  overall  per¬ 
formance  of  the  Hush  House.  Recall  that  the  purpose  of  the  augmenter  tube  is  to  re¬ 
duce  the  exhaust  gas  temperature  and  velocity  to  an  environmentely  acceptable  level 
prior  to  releasing  the  gas  to  the  atmosphere.  The  performance  of  the  Hush  House  is  a 
function  of  the  following  parameters: 

1.  Condition  of  the  flow  field  prior  to  entering  the  engine  inlet 

2.  Magnitude  of  the  mass  flow  rate  of  the  secondary  cooling  air 

3.  Reduction  of  the  gas  temperature  and  velocity  through  the  augmenter  tube 

The  original  augmenter  tube  length  of  21.0  meters  was  reduced  lO'^o  to  a  length  of 
18.9  meters.  The  parametric  study  on  the  effects  of  augmenter  tube  length  on  system 
behavior  used  the  same  k-  c  turbulence  model  as  the  original  model.  If  the  performance 


12 


Figure  3.  Cell  Faces 


of  the  Hush  House  was  not  significantly  reduced  by  a  shorter  augmentcr  tube  then  the 
cost  and  time  of  construction,  as  well  as  the  building  size  could  be  reduced. 

Prior  to  conducting  further  geometrically  based  parametric  studies  the  question 
arose  of  the  elfects  that  dilferent  turbulence  models  had  on  tlic  aerothermal  character¬ 
istics  generated  by  PHOENICS.  'Hie  concern  is  to  know  how  dependent  the  pressure, 
velocity,  and  temperature  fields  were  on  turbulence  modeling.  1  he  GRC.Xl  subroutine 

oilers  the  user  of  PHOENICS  four  turbulence  models:  [Ref  4:  p.  5] 

1.  Algebraic  Model 

2.  Mixing  Length  Model 

3.  k- Mixing  Length  Model 

4.  k-  £  .Model 

Gas  flow  through  the  Hush  House  is  similar  to  How  through  a  three  dimensional 
duct  which  is  a  wall  boundary  problem  between  two  parallel  Y-Z  planes  and  two  parallel 
X-Z  planes.  The  Hush  House  model  is  actually  much  more  complicated  then  just  stated 


Figure  4.  PHOENICS  Flow  Chart 
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Figure  5.  Mush  House  Model  (X-Z  plane) 
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but  in  attempting  to  select  an  existing  turbulence  model  an  assumption  of  the  geometry 
configuration  is  required.  The  algebraic  and  mixing  length  models  as  coded  in  GRLXl 
are  applicable  to  flow  over  a  single  flat  plate  and  therefore  are  not  suitable  for  this 
model.  In  the  k-mixing  length  model  turbulent  viscosity  is  calculated  according  to: 

v,  =  (T-l) 


where  k,  the  turbulence  kinetic  energy  is  a  dependent  variable  solved  for  by  PllOENlCS. 
Q  is  a  constant  set  in  GREXl  and  /„  is  the  mixing  length.  GREXl  has  seven  options 
for  computing  the  mixing  length  l„.  The  first  three  options,  are  applicable  to  flow  over 
a  single  flat  plate  and  are  not  suitable  for  the  Hush  House  model.  The  fourth  option 
selects  the  mixing  length  scale  suitable  for  free  jets.  The  Hush  House  has  a  wall  bounded 
jet.  The  fifth  option  is  suitable  for  free  jets,  but  requires  the  setting  of  nine  constants 
which  can  only  be  determined  through  experimentation  of  the  specific  flow  system.  The 
unavailability  of  these  constants  ruled  this  option  out.  1  he  remaining  two  options  select 
Nikuradse's  mixing  length  scale  which  requires  a  cylindrical  polar  coordinate  grid  system 
for  flow  through  a  cylindrical  duct.  The  k-  c  turbulence  model  is  the  only  PHOENICS 
model  that  is  suitable  for  the  Hush  House  without  drastic  changes  to  grid  geometry. 
However,  since  the  k-  c  turbulence  model  was  used  originally  and  the  other  GREXl 
turbulence  options  are  unsuitable  a  turbulence  model  had  to  be  developed  and  inserted 
into  GROUND.  A  simple  temperature  dependent  turbulent  viscosity  model  was  chosen 
because  of  the  complex  geometry  involed.  The  turbulent  viscosity  term  ENL’T  was  set 
equal  to  a  multiplying  factor  times  the  kinematic  laminar  viscosity  ENUL. 

EMT=  MFx  EML  (3.5) 

ENUL  (commonly  referred  to  as  v,)  is  a  function  of  the  dynamic  viscosity  n  and  the 
density  p  . 
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EXL  L  .  f 


(3.6) 


The  dytiairiic  viscosiiy  of  air  icmpeiaiurc  depcnJcut  over  large  temperature  ranges. 

;i(295'A')  =  18.9  X  10"^  ^  (3.7) 

Ai(2000'a)  =  62.27  x  10~^  ■—  (3.S) 

>71 

A  fourth  order  polynomial  equation  for  (T)  was  developed  by  Andrews  and  Bihlarz  in 
19S1  [Ref  6:  p.C-IJ. 

Ai(7)  =  (.vl  +  .v2  X  74- .r3  X  r  +  x4  x  +  x5  x  T^)  x  \.0  x  10“^  (3.9) 


where: 

xl  =  4.SS56745 
x2  =  5.43232  x  10-^ 
x3  =  -2.4261775  x  10-' 
x4  =  7.9306  X  10-’ 
x5  =  -1.1039S  X  10-'= 

PHOENICS  solves  for  the  temperature  and  density  of  each  cell  therefore  ENUT  is 
temperature-density  dependent  and  is  not  a  constant  value  for  each  cell  in  the  grid  for 
any  given  multiplying  factor, 

Mr) 

EMT^  MFx4^  (3.10) 

The  code  line  ENLT=GRND  in  group  9  of  the  Q1  file  directs  SATELLITE  to  the 
GROL'XD  subroutine  to  evaluate  ENLT  for  each  cell.  The  code  for  evaluating  ENUT 
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is  imbedded  in  group  9.  section  6  of  the  GROL'ND  subroutine  enclosed  in  Appendix 
C. 

The  multiplying  factor  was  arbitrarily  selected  as  200  initially  and  would  be  in¬ 
creased  by  200  after  convergence  of  each  solution.  This  procedure  would  produce  sets 

of  pressure,  velocity,  and  temperature  fields  that  are  based  on  different  kinamatic  tur- 

V5 

bulent  viscosity  values  ranging  from  .0001  to  10  It  was  anticipated  that  for  low 

values  of  kinematic  turbulent  viscosity  the  mixing  effect  of  turbulence  would  be  low  and 
therefore  the  fiow  field  would  be  highly  streamlined.  As  the  turbulent  viscosity  and 
mixing  eflect  increased  it  was  anticipated  that  the  momentum  transport  between  the 
exhaust  gas  and  entrained  cooling  air  would  increase,  resulting  in  lower  velocities  and 
temperatures  exiting  the  e.xhaust  stack. 


IV.  SOLUTION  CONVERGENCE 


A.  SOLITION  METHODS 

PHOENICS  provides  three  methods  for  solving  the  NX  x  NY  x  NZ  number  of 
equations  for  each  dependent  variable: 

1.  Slabwise  Method 

2.  Whole  Field  Method 

3.  Point-By-Point  Method 

A  slab  is  an  array  of  cells  with  the  same  value  of  low-to-high  in  the  Z  coordinate. 
In  the  Hush  House  model  there  are  40  slabs.  The  Slabwise  method  of  solution  solves  for 
all  the  dependent  variables  in  the  lowest  slab  and  then  sequentially  continues  to  the 
highest  slab.  The  sequence  of  solving  equations  from  low  to  high  slabs  is  referred  to  as 
a  sweep.  The  solution  of  the  discretized  equation  for  any  dependent  variable  requires 
values  from  each  neighboring  cell.  In  the  Slabwise  solution  the  equations  for  values  in 
a  slab  require  values  from  the  slabs  below  and  above  it.  Adjustments  to  the  values  in 
the  ne.xt  higher  slab  will  invalidate  the  values  just  solved  for  in  the  lower  slab,  therefore 
many  sweeps  are  required  until  all  the  discretized  equations  balance  and  no  further  ad¬ 
justments  to  neighboring  cells  are  required.  (Ref  7;  p.  2.5] 

The  Whole  Field  Method  of  solution  simultaneously  solves  the  NX  x  NY  x  NZ 
number  of  equations  for  each  dependent  variable  in  one  sweep  thereby  reducing  the 
number  of  sweeps  compared  with  the  Slabwise  Method  by  a  factor  of  NZ.  The  Whole 
Field  .Method  requires  much  more  computer  storage  space  then  the  Slabwise  Method 
and  should  not  be  used  if  storage  space  is  limited.  [Ref.  7;  p.  2.5] 

The  Point-By-Point  Method  of  solution  does  not  employ  PHOENICS  simultaneous 
linear  equation  solver  for  the  evaluation  of  the  dependent  variables.  The  node  values 


of  each  variable  arc  updated  by  simple  arithmetic  substitution.  This  simple  substitution 
reduces  the  amount  of  computational  work  but  requires  many  cycles  to  reach  a  final 
solution. 

B.  METHODS  OF  DL  FERMINING  CONVERGENCE 

.4  converged  solution  is  one  in  which  the  imbalance  between  the  left  and  right  hand 
sides  of  every  discretized  equation  is  reduced  to  an  acceptable  small  value.  Due  to  the 
iterative  procedure  utilized  by  PHOENTCS  in  solving  the  discretized  equations,  conver¬ 
gence  is  not  guaranteed.  User  intervention  is  usually  required  to  obtain  a  converged 
solution. 

There  are  four  means  available  to  the  user  of  PHOENTCS  to  evaluate  whether  a 
solution  is  converging  or  diverging: 

1.  Tabulation  of  residuals  vs  sweep 

2.  Conservation  of  mass  flow 

3.  Conservation  of  energy 

4.  Comparison  of  successive  sweeps 

Residuals  calculated  for  each  dependent  variable  are  the  sum  of  the  absolute  values  of 
the  discretized  equation  imbalance  from  each  ceil.  The  decrease  in  residuals  with  each 
successive  sweep  is  the  primary  means  of  determining  convergence.  The  output  data  file 
generated  by  PHOENTCS  provides  the  values  of  residuals  vs  sweeps  for  each  dependent 
variable  in  tabular  and  graphical  form.  Figure  9  on  page  28  show  graphical  examples 
indicating  monotonic  convergence.  Monotonic  convergence  is  characterized  by  sequen¬ 
tial  decrease  of  residuals  vs  sweeps.  Monotonic  convergence  will  usually,  but  not  al¬ 
ways,  yield  a  final  solution  in  the  shortest  amount  of  time.  Monotonic  convergence  is 
not  always  possible  or  may  be  very  difficult  to  obtain.  Residuals  may  have  a  damped 
oscillatory  decay  as  shown  in  Figure  10  on  page  29.  The  computational  time  required 
to  reach  a  converged  solution  with  the  residuals  converging  this  way  may  be 
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considerably  greater  than  if  nionotonic  convergence  could  be  obtained.  Ho\\c\er  a 
rapidly  decreasing  damped  oscillatory  decay  can  converge  faster  than  o  slowlv  decreasing 
monotonic  decay.  Divergence  of  a  solution  is  determined  by  the  residuals  increasing 
with  the  number  of  sweeps  as  shown  in  Figure  11  on  page  30.  A  diverging  solution  if 
unaltered  will  most  likely  result  in  "self  destruction"  of  the  simulation. 

There  are  no  cut  and  dry  methods  for  obtaining  and  maintaining  a  converged  sol¬ 
ution.  The  general  guidelines  for  obtaining  a  converged  solution  are  discussed  in  part  2 
section  S  of  Reference  7.  In  choosing  a  viable  set  of  convergence  control  options  from 
the  large  number  of  possibilities,  the  PHOENICS  programmer  must  use  contmon  sense, 
experience,  and  trial  and  error. 

Three  convergence  control  options  were  utilized  in  the  Hush  House  model: 

1.  Under  relaxation  factor 

2.  Termination  criteria  for  inner  iterations 

3.  Incrementation  of  source  terms 

The  under  relaxation  factor  controls  the  rate  of  change  of  a  dependent  variable  from 
sweep  to  sweep.  Due  to  the  iterative  procedure  utilized  by  PHOENICS,  improper  set¬ 
ting  of  the  under  relaxation  factor  of  a  dependent  variable  can  result  in  a  diverging  sol¬ 
ution.  Under  relaxation  factors  range  from  just  greater  than  zero  to  one.  Small  values 
of  the  under  relaxation  factor  require  more  iterations  for  convergence  to  be  obtained  and 
thus  more  computer  time  to  reach  a  final  solution.  An  optimum  under  relaxation  factor 
is  one  that  minimizes  computer  time  and  also  prevents  divergence.  The  determination 
of  an  optimum  under  relaxation  factor,  as  in  the  choice  of  convergence  control  options, 
requires  experience  and  trial  and  error.  The  under  relaxation  factor  is  activated  by  the 
use  of  the  RELAX  command  in  group  17  of  the  Q1  file. 

The  second  convergence  option,  termination  criteria  for  inner  iterations,  controls  the 
number  of  iterations  performed  on  a  dependent  variable  for  each  sweep.  A  larger 
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number  of  iterations  per  sweep  will  yield  a  more  accurate  value  for  that  variable  and  thus 
result  in  fewer  sweeps  but  the  computing  time  per  sweep  is  increased.  A  balance 
between  the  number  of  inner  iterations  per  sweep  and  the  number  of  sweeps  is  required. 
It  would  be  a  waste  of  computing  time  to  iterate  on  a  variable  to  a  degree  of  accuracy 
exceeding  the  required  accuracy  of  the  output  The  number  of  inner  iterations  is 
controled  by  the  command  LITER  in  group  16  of  the  Q1  file. 

The  incrementation  of  source  terms  is  a  method  used  to  slowly  apply  the  excitation 
to  the  flow  system  thereby  reducing  the  jump  required  in  solving  the  conservation 
equations.  The  total  excitation  to  the  system  is  obtained  by  applying  the  source,  such 
as  mass  flow  rate  or  temperature,  in  small  increments.  The  magnitude  of  the  excitation 
is  increased  after  run  of  sweeps  has  converged  based  on  the  previous  excitation  magni¬ 
tude.  This  procedure  is  in  effect  a  transient  to  steady  state  How  condition  solution. 

For  the  first  parametric  study  involving  the  effects  of  augmenter  tube  length  on 
system  behavior,  a  converged  solution  with  a  ten  percent  reduction  in  augmenter  tube 
length  was  desired.  Having  little  e.xperience  at  this  time  with  PHOENICS.  the  same 
under  relaxation  factors  and  number  of  inner  iterations  were  used  as  in  the  original  Hush 
House  model  Q1  file  developed  by  Nicolaus.  The  variable  FACT  is  used  to  define  the 
magnitude  of  the  FALSDT  under  relaxation  factor.  A  value  of  I'ACT  =  0.75  was  ap¬ 
plied  to  the  dependent  variables  L  I,  VI,  W'l,  KE,  and  EP.  A  value  of  2.5  x  FACT  was 
used  for  the  dependent  variable  HI.  The  variable  LINRLX  is  used  to  define  the  linear 
under  relaxation  factor.  A  value  of  LINRLX  =  0.3  was  used  for  the  variable  PI.  The 
default  value  of  LITER  equal  to  20  iterations  per  sweep  was  used  for  all  dependent 
variables.  The  DF09  file  from  the  original  Hush  House  converged  solution  was  re¬ 
started.  To  avoid  a  jump  or  shock  to  the  system  the  original  augmenter  tube  length  was 
initially  reduced  by  5  percent  with  2500  sweeps  required  for  a  converged  solution.  The 
augmenter  tube  length  was  shortened  an  additional  5  percent  of  the  original  length  with 
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3500  sweeps  required  for  a  converged  solution.  Each  of  the  32  runs,  consisting  of  either 
300  or  SO  sweeps,  exhibited  monotonic  convergence. 

For  the  second  parametric  study  involving  the  effects  of  turbulent  viscosity  on  sys¬ 
tem  behavior  a  trial  and  error  procedure  for  obtaining  convergence  was  required.  The 
multiplying  factor  .\1F  of  equation  (3.10)  was  set  equal  to  200.  The  DF09  file  from  the 
original  Hush  House  model  was  used  for  a  restart  with  the  exceptions  of  the  dependent 
variables  KE  and  EP  which  were  not  required  for  the  new  temperature  dependent  tur¬ 
bulence  model  .  Still  having  limited  experience  with  PHOENTCS,  the  same  under  relax¬ 
ation  factors  and  number  of  inner  iterations  were  used  as  in  the  original  Hush  House 
model  Q1  file.  The  above  procedure  resulted  in  divergence  of  the  solution.  The  number 
of  iterations.  LITER,  for  the  pressure  variable  PI  was  increased  from  the  default  value 
of  20  to  40  because  in  the  SI.MPLE  algorithim  the  pressure  is  solved  for  first  from  the 
continuity  equation  and  then  the  velocity  components  are  solved  from  the  momentum 
equations  which  have  a  pressure  gradient  term.  A  more  accurate  value  of  pressure  ob¬ 
tained  by  more  iterations  will  yield  more  accurate  velocity  components  without  ha\ing 
to  increase  the  number  of  iterations  performed  on  them.  The  variable  LINRLX  was 
reduced  from  0.3  to  O.l  for  pressure  variable  PI.  The  under  relaxation  factor.  F.ACT 
was  reduced  from  0.75  to  0.25  for  the  remaining  variables  Ul.  \T.  \V1,  HI.  These  three 
changes  resulted  in  an  oscillation  of  residuals  with  an  overall  mean  effect  of  neither  in¬ 
creasing  nor  decreasing.  It  was  evident  that  after  nearly  a  thousand  sweeps  that  the 
solution  was  not  going  to  converge  within  an  acceptable  time  frame.  The  next  attempt 
at  obtaining  a  converged  solution  involved  incrementing  the  temperature  source  from 
the  jet  engine  exhaust.  The  exhaust  temp  was  set  equal  to  ambient  air  at  295  °  K  and 
the  problem  was  started  from  scratch  and  not  a  restart  from  the  original  model.  The 
mass  flow  rate  was  maintained  at  3S.49  kg  s.  The  first  few  runs,  consisting  of  80  sweeps 
each,  exhibited  monotonic  convergence.  FACT  was  then  increased  from  0.25  to  1.0  to 


minimize  computer  CPU  time.  Monotonic  convergence  was  present  during  t!ie  16S() 
sweeps  required  for  a  final  solution  with  engine  exhaust  temperature  at  295  °  K.  Engine 
exhaust  temperature  was  increased  to  500  “  K  and  a  converged  solution  was  obtained  in 
52S0  sweeps.  With  the  engine  exhaust  temperature  incremented  to  700  °  K.  an  addi¬ 
tional  3100  sweeps  were  required  for  a  converged  solution.  With  the  engine  exhau'-t 
temperature  set  to  900  °  K,  the  residuals  oscillated  severly,  yet  the  mean  values  were 
neither  converging  nor  diverging  over  successive  runs.  FACT  was  incremently  reduced 
to  0.1  and  LINRLX  remained  at  0.1.  After  nearly  5000  sweeps,  oscillation  of  the  resi¬ 
duals  continued  with  the  same  mean  effect  as  exhibited  by  the  residuals  of  the  first  few 
runs.  This  dilemma  of  non  convergence  non  divergence  is  a  prime  example  of  not 
choosing  the  correct  convergence  control  options.  It  should  be  fully  understood  that  the 
solution  of  the  flow  field  with  engine  exhaust  temperature  at  900  °  K  is  feasible  but  it 
requires  a  different  approach  than  that  used  for  the  lower  temperatures.  Due  to  the 
shortage  of  time  remaining  for  this  turbulence  study  and  the  anticipated  convergence 
problems  at  higher  temperatures,  it  was  decided  to  conduct  the  parametric  study  of  tur¬ 
bulence  using  the  converged  solution  with  engine  exhaust  temperature  at  700  °  K. 

Three  additional  turbulent  viscosity  models  were  simulated  with  values  of  the  mul¬ 
tiply  factor  .\1F  in  equation  (3.10’/  set  to  lOuO,  20000,  and  200,000.  In  addition,  a  sol¬ 
ution  for  the  original  Hush  House  model  (k-  c  turbulence  model)  with  the  jet  engine 
exhaust  temperature  set  700  °  K  was  simulated. 


27 


RESIDL'ALSrNORMALIZED) 


RESIDUALS  vs  SWEEPS 


0.50 


0.10 


HU 
U  U 
H  U 
H 

M  U 
H 

H  H 

H 

V  H 

V  M  HH 
.P  U  H 

H 

V  M  HH 

U  HH 

W  H 

V  H  H 

P  U  HH 

W  HH 

V  U  HH 

M  H 

U  HH 

P  V  W  HH 

V  U  HH 

V  W  U  H  H 

VV  W  U  HH 

P  V  V  UU  H 

HU  HH 

P  WV  UUU  ri  H 

P  WW  UU  HH 

P  H  UU  U  HH 

P  H  H  UU  HH 

P  WW  UU  HH 

VWW  UUU  HH 

P  VWW  UUU  H 

p  V  WW  UUUU 

P  WWV  UUU 

P  WWW  UUU 

P  WW  WW  U  UU  HH 

PP  MW  WV  UUU  HH 

PP  WWW  VWV  UUU  U  HH 

PP  WW  W  VW  UUU  HHH 

PP  WWW  VWW  UUU  HH 

PPP  P  WWW  W  VW  V  UUUU  HH 

PPP  WWWWW  WWV  UUU  HH 

PPPP  PP  WW  WWWWWWV  VWW  UUUU  HH  H 

PPPPP  PPPP  W  WWWWWWV  WV  UUUU  HH 

PPP  PPPPPPPP  PP  H  WUWWWWWW  UUUHMH 

PPPPP  PPPPPPP  PPPP  PP  WWWWWWH  H 

. • . • . • . • . • . * . • . •  .  .  .P.  .P  .PPPPPPW.WKHH 

.1  .2  .5  .5  .7  .8  .9  1.0 


HH 


HHH 


HH 


SWEEPS(NORMALIZED) 


Figuie  9. 


Sample  Output  Residuals  vs  S^veeps  (Converging  Solution) 
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Figure  10.  Sninple  Output  Residuals  vs  S^eep  (Converging  Solution) 
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Figure  1 1.  Sample  Output  Residuals  vs  Sueeps  (Diverging  Solution) 
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V.  USING  THE  PHOENICS  CODE 

The  efTicient  and  successful  utilization  of  PHOENICS  requires  a  thorough  under¬ 
standing  of  the  PHOENICS  Code,  the  user  manual,  and  the  various  computing  facilities 
required  to  operate  the  PHOENICS  Code.  It  is  extremely  difficult  to  obtain  a  thorough 
understanding  of  PHOENICS  in  a  six  month  period.  User  manuals  in  general  can  not 
provide  every  bit  of  information  required  to  fully  operate  and  understand  a  system.  This 
section  is  provided  as  an  additional  resource  of  information  to  help  eliminate  conunon 
problems  encountered  using  PHOENICS.  It  is  intended  that  this  information  will  re¬ 
duce  the  many  setbacks  that  the  beginning  users  of  PHOENICS  arc  sure  to  stumble  on. 

A.  ALLOCATING  MVS  DATA  SETS 

The  PHOENICS  code  is  run  as  a  batch  job  on  the  IBM  3033  main  frame  computer. 
The  batch  system  at  the  NPS  Computing  Center  is  MVS  (Multiple  Virtual  System).  The 
MVS  system  introduces  programs  to  the  computer,  initiates  their  execution,  and  sched¬ 
ules  all  resources  and  services  required.  MVS  data  sets  must  be  allocated  prior  to  run¬ 
ning  the  Q1  file.  For  cartesian  coordinate  system  grid,  the  following  data  sets  are 
required: 

1.  MSS.S'”--’''.PHOENICS.DF09 

2.  MSS.S'-'-'.PHOENTCS.DElO 

3.  MSS.S-*--.PHOENTCS.LOAD 

4.  MSS.S-'-NPHOENICS.SRC 

Where  should  be  replaced  with  the  users  NPS  computer  ID  number.  .MVS  data 
sets  are  allocated  by  use  of  the  CMS  command  .MVSHELP.  Select  option  2  for  allo¬ 
cation  of  new  data  sets.  The  computer  will  prompt  for  the  required  information.  If  the 
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body  fitted  coordinate  (BFC)  grid  system  is  used  instead  of  the  cartesian  coordinate  grid 

system,  then  two  additional  data  sets  are  required: 

1.  MSS.S'--'\PMOEMCS.DF12 

2.  MSS.S----.PHOFNICS.DF16 

Table  1  lists  the  data  set  parameters  for  the  six  required  files. 


Table  1.  DATA  SET  PARAMETERS 


DATA 

SET 

RECFM 

LRECL 

BLKSl 

DSORG 

DF()9 

VBS 

19004 

19004 

sequential 

DFIO 

F 

Soft 

SOO 

sequential 

DF12 

F 

400 

400 

sequential 

DFl  6 

F 

4(»0 

400 

sequential 

LOAD 

U 

N  A 

13030 

partitioned 

SRC 

FB 

SO 

so 

partitioned 

Another  option  of  allocating  .\1VS  data  sets  and  by  far  the  easiest  method  is  to  re¬ 
name  a  previous  PHOENICS  users  data  sets  using  the  CMS  command  MVSHELP  op¬ 
tion  4  the  rename  option.  All  MVS  data  sets  are  accessible  to  anyone  via  the  C.MS 
commands  MVSHELP  and  GETMVS.  This  procedure  will  eliminate  the  problem  of 
allocating  space  for  empty  data  sets  for  PHOENICS. DFl 2  and  PHOENICS. LO.-\D 
which  results  in  error  messages  when  running  the  Q1  file  for  the  first  time. 

A  few  words  of  caution  are  needed  at  this  point  concerning  MVS  data  sets.  If  tlie 
DF09  file  is  transferred  from  another  user,  the  restart  option  in  group  24  of  the  Q1  file 
should  not  be  used  (i.e.,RESTRT  should  be  deactivated)  the  first  time  the  Q1  file  is  run 
unless  the  data  stored  on  the  DF09  file  is  to  be  used  as  input  data.  After  the  first  run 
invoke  the  restart  option.  DF09  is  the  file  that  the  grid  geometry  and  PHOENICS 
generated  output  is  stored  on.  Each  time  the  problem  is  run  the  output  from  the  previ¬ 
ous  run  is  read  from  the  DF09  file  as  the  input  data  for  the  present  run.  The  output  for 
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any  run  writes  over  the  previous  output  on  the  DF09  file.  file.  This  procedure  allows 
sequential  runs  to  be  performed  with  the  latest  output  being  submitted  as  the  present 
input  data  for  a  restart.  In  order  to  save  good  output  files  and  prevent  them  from  being 
written  over,  the  DF09  file  must  be  copied  onto  a  different  backup  MVS  data  set.  The 
program  COPYTAPE  listed  m  Appendix  B  provides  the  means  of  copying  any  MVS 
data  set  to  another  MVS  data  set.  The  new  MVS  Data  File  must  be  allocated  prior  to 
copying  and  have  the  same  MVS  parameters  as  the  one  that  is  being  copied.  The  name 

of  the  new  file  is  arbitrary  but  the  recommended  name  form  is: 
MSS.S-*^-.PHOENTCS.BACKLP 

Where  is  the  users  computer  ID  number  and  BACKUP  is  the  new  file  name  (or 
any  other  file  name).  The  function  of  COPYTAPE  is  to  copy  the  SYSUTl  data  set  to 
the  SYSUT2  data  set.  COPYTAPE  can  be  used  to  copy  the  PHOENICS  output  file 
onto  the  backup  file  or  copy  the  backup  file  onto  the  DF09  file  which  would  be  used  as 
the  input  data  for  the  ne.xt  run  of  the  Q1  file.  It  is  strongly  recommended  to  copy  the 
DF09  file  onto  a  backup  file  on  regular  intervals  such  as  on  a  daily  or  bi-weekly  basis. 
It  is  possible  that  the  DF09  file  could  contain  bad  data  as  a  result  of  a  change  in  the 
Q1  file  or  the  submission  of  a  dificrent  Q1  file  having  the  same  computer  ID  numbers. 
If  this  dilemma  occurs  and  the  backup  files  are  utilized,  only  a  loss  of  several  days  of 
computing  would  be  lost,  vice  weeks  or  months  if  backup  files  were  not  utilized. 

B.  IMPLEMENTING  GROUND 

GROUND  is  the  program  for  interaction  with  the  EARTH  solution  program.  The 
implementation  of  GROUND  requires  five  steps: 

1.  Insert  relevent  equations  into  the  file  GROUND  FORTIU^N 

2.  Subi-nit  COPYGRD 

3.  Submit  COMPGRD 

4.  Compact  LOAD  and  SRC  MVS  data  sets 
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5.  Submit  LINKEAR  JCL 


Group  9  of  the  GROUND  FORTR.-W  file  is  shown  in  Appendix  C.  Section  5  contains 
the  coding  required  to  evaluate  ENLT  (equation  (3.10))  used  in  the  Mush  House  tur¬ 
bulence  study.  The  programs  COPYGRD,  COMPGRD,  AND  LINKE.AR  are  included 
in  Appendix  D  for  reference.  Compacting  the  LOAD  and  SRC  files  requires  the  se¬ 
quential  submission  of  three  programs  for  each  file.  For  the  LD.^D  file  the  following 
programs  listed  in  Appendix  E  are  required: 

PHOENXOl  LOADLIB 
PHOENX02  LOADLIB 
PHOENX03  LOADLIB 

For  the  SRC  file  the  following  programs  listed  in  Appendix  E  are  required: 

PHOENXOl  SOURCELB 
PHOENX02  SOURCELB 
PHOENX03  SOURCELB 

When  the  LOAD  data  set  has  been  copied  from  a  previous  user  it  will  contain  the 
equations  placed  in  the  GROUND  FORTR.AN  file  by  the  previous  user.  If  the 
equations  in  GROUND  FORTR^AN  are  not  required,  ensure  that  the  command 
L'SEGRD=  F  is  inserted  in  group  19  of  the  Q1  file. 

C.  DATA  TRANSFER  FROM  MVS  TO  VAX 

Although  the  PHOENTCS  Code  is  available  on  the  NPS  computer  center  IBM 
Mainframe,  the  graphics  program  PHOTON  is  only  available  on  the  VAX  system. 
Therefore  in  order  to  have  data  transfer  of  the  DF09  file  from  MVS  to  the  VAX  2000 
system  a  data  transfer  operation  is  required.  Presently  there  are  two  methods  available 
for  data  transfer  at  NPS.  The  first  method  involves  copying  the  DF09  data  set  to  a  CMS 
file  and  then  transferring  the  C.MS  file  via  a  modem  from  the  mainframe  to  the  VAX 
system.  The  computer  hardware  and  software  for  the  procedure  has  not  been  fully  de¬ 
bugged.  More  often  than  not,  the  transfer  of  data  is  interrupted  and  data  transfer  is 
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seldom  successful.  Until  the  system  has  been  debugged,  it  is  strongly  recommended  to 
avoid  this  procedure  as  it  is  simply  inefficient  and  undependable. 

The  second  method  available  is  to  transfer  the  DF09  MVS  data  set  to  a  9  track  tape 
and  then  transfer  from  the  tape  to  a  VAX  cartridge.  Prior  to  transferring  the  DF09  data 
file  to  the  9  track  tape,  the  unformatted  DF09  file  must  be  formatted  so  that  the  file  can 
be  read  by  PHOTON.  The  program  FORM09  in  Appendix  F  transforms  the  original 
unformatted  DF09  file  into  a  new  formatted  DF09C  file.  The  new  DF09C  file  must  be 
allocated  using  MVSHELP  prior  to  submitting  FORM09.  The  proper  setting  of  the 
parameters  for  DF09  are  crucial.  The  data  set  organization,  record  format,  logical  re¬ 
cord  length,  and  blocksize  must  be  set  as  listed  below.  The  number  of  primary  and 
secondary  storage  tracks  allocated  to  the  DF09C  file  may  vary  depending  on  the  size  of 
the  DF09  file. 

Data  set  organization  =  Sequential  (PS) 

Record  Format  =  Fixed  Block  (FB) 

Logical  record  length  =  80 
Blocksize  =  800 

Primary  Allocation  =  250  (TRK) 

Secondary  =  16  (TRK) 

A  nine  track  tape  may  be  checked  out  from  the  NTS  computer  center.  A  non  la¬ 
beled  tape  with  a  density  of  either  1600  or  6250  bits  per  inch  (bpi)  is  required.  If  several 
files  are  to  be  transfered  to  tape,  the  6250  density  tape  is  recommended  due  to  its  larger 
capacity.  The  program  SEQ5,  a  utility  program  available  through  the  MVSHELP 
command  ,  copies  the  formatted  DF09C  file  to  the  9  track  tape.  The  tape  to  VAX  car¬ 
tridge  transfer  equipment  is  located  in  the  lab  of  the  NTS  Meteorology  Department 
(Donna  Bur\  ch  is  the  point  of  contact  for  computer  access  authorization).  The  follow¬ 
ing  VAX  commands  arc  required  for  data  transfer  from  tape  to  VVXX  cartridge; 
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Mount  For  Record  =  80  Blcoksize  =  800  MS  AO: 

Mount  Record  =  80  Blocksize  =  800  MUAO: 

Copy  MSAO:  MUAO:  FX.FT 
Dismount  MSAO: 

Dismount  MUAO; 

During  the  data  transfer  ignore  the  "incompatible  attributes"  message.  It  has  no  eflect 
on  the  transfer  operation.  When  using  the  cartesian  coordinate  grid  system  (  which  does 
not  reqiure  the  use  of  .MVS  DF12  data  set  ),  the  DF09  file  on  the  VAX  cartridge  does 
not  require  any  editing  to  be  compatible  with  the  Phoenics  version  1.4  on  the  VAX. 

In  running  the  PHOENICS  Code  it  is  inevitable  that  computing  errors  will  occur. 
Listed  below  are  some  of  the  common  error  codes  encountered  along  with  the  corrective 
action  required: 

Abend  code  322  -  allocated  CPU  time  has  been  exceeded.  The  number  of  sweeps 
per  run  must  be  reduced.  In  the  Hush  House  model  a  maximum  of  90  sweeps  for  a  15 
ntinute  class  G  job  and  360  sweeps  for  a  60  minute  class  J  job  were  possible.  The  time 
per  sweep  is  dependent  on  the  number  of  equations  being  solved  and  the  number  of  it¬ 
erations  (LITER)  for  each  variable. 

Abend  Code  SE37  -  The  LOAD  and  SRC  files  require  compacting.  The  compacting 
procedure  is  discussed  in  section  B  of  this  chapter. 


VI.  COMPUTER  MODEL  RESULTS 


The  parametric  study  of  the  aerothermodynamics  of  the  NAS  Jackson\ille,  Florida 
Hush  House  consisted  of  two  independent  studies.  The  first  study  investigated  the  ef¬ 
fects  of  a  10'!  b  reduction  of  the  augmenter  tube  length  on  system  behavior.  The  second 
study  investigated  the  elFects  of  turbulence  model  on  system  behavior.  For  ease  in  dif¬ 
ferentiating  the  diflerent  Hush  House  models,  the  following  terms  are  defined  and  will 
be  referred  in  the  discussion  of  results: 


Model  1  -  Original  Hush  House  Model  utilizing  the  k  -  £  turbulence  model  (TJET  = 
1944  ^  K) 

Model  2  -  Hush  House  Model  with  lOlb  reduced  augmenter  tube  length  utilizing  the 
k  -  e  turbulence  model  {TJET  =  1944  ^  K) 

Model  3  -  Hush  Flouse  Model  utilizing  equation  (3.10)  for  turbulent  viscositv  (ME 
=  200000,  TJET  =  700  K)  “ 

Model  4  -  Hush  House  .Model  utilizing  equation  (3.10)  for  turbulent  viscositv  (MF 
=  2U000.  TJET  =  700  K) 

Model  5  -  Hush  House  Model  utilizing  equation  (3.10)  for  turbulent  viscositv  (ME 
=  1000,  TJET  -  700  K) 

Model  6  -  Hush  House  Model  utilizing  equation  (3.10)  for  turbulent  viscositv  (MF 
=  200.  TJET  =  700  K)  ^ 

Model  7  -  Original  Hush  House  Model  utilizing  k  -  e  turbulence  model  (TJET  -  700 
K) 

The  tabulated  numerical  results  generated  by  PHOENICS  for  each  of  the  7  models 
are  comprised  of  nearly  10000  lines  of  data.  For  ease  of  comparing  results,  output  data 
is  presented  in  graphical  form  for  pressure,  temperature,  and  velocity  profiles  and  con¬ 
tours. 

A.  AUGMENTER  Tl^BE  LENGTH  STUDY 

This  study  compares  the  results  of  model  1  (the  original  model)  and  model  2  (the 
original  model  with  a  shortened  augmenter  tube).  The  pressure,  temperature,  and  ve- 
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locity  profiles  from  IZ=  1-40  at  IX=  1  and  IY  =  4  for  model  1  and  model  2  are  shown 
in  Figure  12  on  page  43,  Figure  13  on  page  44  and  Figure  IJ  on  page  45.  The  pressure 
profile  shows  identical  values  of  pressure  for  both  models  from  Z  =  0  to  Z=  27.4  meters. 
The  pressure  in  model  2  is  approximately  40  pascals  higher  than  in  model  1  from 
Z=27.4  to  Z  =  45.6  meters.  The  latter  distance  represents  the  107  b  shorter  augmenter 
tube  length.  The  remaining  section  of  the  model  2  profile  is  identical  to  the  model  1 
profile  except  that  it  is  shifted  to  the  left  a  distance  of  2.1  meters,  which  is  equal  to  the 
reduction  in  augmenter  tube  length.  The  distance  axis  is  normalized  by  using  the  cell 
numbers  in  the  Z  direction  vice  actual  distance  in  the  Z  direction  as  shown  in 
Figure  15  on  page  46.  The  normalized  profile  shows  more  clearly  that  there  is  little  if 
any  change  in  pressure  magnitudes  between  the  two  models.  Similarly  for  the  normal¬ 
ized  temperature  and  velocity  profiles  shown  in  Figure  16  on  page  47  and  Figure  17  on 
page  48  there  is  no  significant  changes  in  temperature  and  velocity  magnitudes  between 
the  two  models.  Velocity  vector  plots,  temperature  contours,  and  pressure  contours  of 
model  1  and  model  2  are  shown  in  Figure  18  on  page  49  through  Figure  29  on  page 
60.  Figure  18  on  page  49  and  Figure  19  on  page  50  show  the  velocity  vector  plots  in  the 
test  bay  at  engine  level  (IY  =  4)  are  identical  for  model  1  and  model  2.  Figure  20  on  page 
51  and  Figure  21  on  page  52  show  the  vector  velocity  plots  at  the  exhaust  of  the 
augmenter  tube  are  identical  for  model  1  and  model  2.  Figure  22  on  page  53  and 
Figure  23  on  page  54  show  the  temperature  contours  in  the  entrance  of  the  augmenter 
tube  are  identical  for  model  1  and  model  2.  Figure  24  on  page  55  and  Figure  25  on  page 
56  show  the  temperature  contours  at  the  exhaust  of  the  augmenter  tube  arc  identical  for 
model  1  and  model  2.  Figure  26  on  page  57  and  Figure  27  on  page  58  show  the  pressure 
contours  at  the  entrance  of  the  augmenter  tube  are  identical  for  model  1  and  model  2. 
Figure  28  on  page  59  and  Figure  29  on  page  60  show  the  pressure  contours  at  the  ex¬ 
haust  of  the  augmenter  tube  are  identical  for  model  1  and  model  2.  These  figures  pro- 
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vide  verificalion  that  the  shorter  augnienter  tube  has  little  efTect  on  the 
aerothermodynamics  of  the  How  field.  In  addition  to  comparing  profiles  and  contours 
of  model  1  and  model  2,  the  augmentation  ratios  of  each  model  are  compared.  The 
augmentation  ratio  is  defined  as  follows: 


A.R.  = 


■^^/  V  ~  ^^E\0 

^hsG 


where: 


M.v  =  inlet  mass  flow  rate 


^^ENG  “  engine  mass  flow  rate 


The  inlet  mass  flow  rate  is  the  sum  of  the  mass  flow  rates  through  the  five  air  inlet 

Ac 

doors  in  the  test  bay.  The  engine  mass  flow  rate  of  3S.49  is  the  same  in  both 
models.  The  augmentation  ratio  is  a  measure  of  the  efficiency  of  the  Hush  House  in 
regards  to  the  magnitude  of  the  mass  flow  rate  of  secondary  cooling  air  that  can  be 
drawn  into  the  Hush  House  by  the  "entrainment  elTect"  of  the  jet  engine  exhaust.  The 
higher  the  augmentation  ratio  the  higher  the  efficiency  of  the  Hush  House  design.  The 
augmentation  ratio  for  model  1  is  10.66  and  for  model  2  is  10.87.  This  amounts  to  an 
increase  of  2%  of  the  augmentation  ratio  for  model  2  compared  to  model  1. 

The  10?'o  shorter  augmenter  tube  is  slightly  more  efficient  than  the  original  model. 
Although  the  pressure,  temperature,  and  velocity  of  the  exhaust  gas  as  it  enters  the  at¬ 
mosphere  is  not  effected  by  the  shorter  tube,  the  noise  level  (not  analysed  in  this  thesis) 
is  expected  to  be  higher.  Recall  that  the  augmenter  tube  is  lined  with  acoustic  absorp¬ 
tion  material.  Reducing  the  augmenter  tube  reduces  the  surface  area  of  the  a  :>sorption 
material  therebv  reducing  the  amount  of  acoustic  enerev  absorbed. 
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B.  TURBULENCE  STUDY 

The  elTects  of  turbulence  model  on  system  aerothermal  behavior  are  presented  in 
pressure,  temperature  and  velocity  profiles  and  contours.  The  magnitudes  of  the  turbu¬ 
lent  viscosity,  resulting  from  equation  (3.10)  ,  in  the  arbitrarily  selected  cell  defined  by 
IX  =  1,  IY  =  4,  and  IZ=20  for  model  3,  model  4,  model  5,  model  6.  and  model  7  are  as 
follows: 


model  3 

MF=  200.000 

vr^6.4^ 

model  4 

MF=  20.000 

.^1  I 

vv~l.l  — r 
nr 

model  5 

MF=  1000 

vv:t0.06^ 

model  6 

.MF=200 

vv^tO.Ol 

//j- 

model  7 

k-  £  model 

Vj.~4.4  -—7 

It  should  be  stressed  that  these  values  of  iv  are  representative  values  of  turbulent 
viscosity  in  one  cell  only.  Values  of  Vj-  are  different  in  each  of  the  cells  of  the  model. 
The  results  from  model  5  and  model  6  are  essentially  the  same  and  are  therefore  pre¬ 
sented  as  a  single  model  referred  to  as  model  5&6.  Model  7  (k-  c  model)  is  included  for 
the  purpose  of  providing  a  relative  reference  point  to  gage  the  effects  of  the  temperature 
dependent  turbulence  models.  The  pressure  profiles  for  model  3.  model  4.  model  5&6, 
and  model  7  are  shown  in  Figure  33  on  page  64.  In  the  test  bay  area  prior  to  the  engine 
inlet  (Z  distance  =  0  -  17.6  meters)  the  pressure  profiles  for  each  model  are  appro.xi- 
mately  constant.  The  noticable  trend  is  that  as  the  turbulent  viscosity  is  increased  from 
vv~0.01  to  vv5:6.4  the  constant  pressure  profile  for  each  model  decreases  by  roughly  20 
pascal.  The  trend  continues  for  the  large  pressure  depression  just  prior  to  the  engine  inlet 
(Z  distance  =  18.6  meters).  After  the  engine  exhaust  (Z  distance  =  22.3  meters),  little 
correlation  between  the  pressure  profiles  and  the  magnitude  of  turbulent  viscosity  is 
observable.  For  low  magnitudes  of  turbulent  viscosity  (17:^0.06  and  iv^tl.l  )  the  pres¬ 
sure  profile  through  the  augmenter  tube  shows  little  variation.  For  higher  magnitudes 
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of  turbulent  viscosity  (v>5:4.4  and  \v2:6.4  )  the  pressure  profile  through  the  augmenter 
tube  varies  from  a  nagative  pressure  to  a  positive  pressure. 

The  temperature  profiles  for  model  3.  model  4,  model  5&6.  and  model  7  are  shown 
in  Figure  30  on  page  61.  The  temperature  profiles  of  the  four  models  are  identical  in 
the  test  bay  (295  °  K),  through  the  engine  (250  °  K)  to  the  engine  outlet  (700  °  K).  It 
is  clearly  shown  that  as  the  magnitude  of  turbulent  viscosity  increases,  the  temperature 
decrease  through  the  augmenter  tube  increases  as  a  result  of  the  increased  mixing  effect. 

The  velocity  profiles  for  model  3,  model  4,  and  model  5&6,  and  model  7  are  shown 
in  Figure  34  on  page  65.  The  velocity  profiles  are  similar  to  the  temperature  profiles  in 
that  the  velocity  decrease  through  the  augmenter  tube  is  increases  as  the  magnitude  of 
the  turbulent  viscosity  increases. 

Velocity  vector  plots  in  the  test  bay  for  model  3,  model  4,  and  model  5&6  are  shown 
in  Figure  35  on  page  66.  Figure  36  on  page  67,  and  Figure  37  on  page  68.  The  velocity 
vector  plot  for  model  7  is  not  included  in  the  figures  becuase  of  its  similarity  with  the 
vector  plot  of  model  4.  The  velocity  vectors  in  model  7  have  the  same  direction  but  a 
magnitude  1.5  times  greater  than  the  velocity  vector  plots  in  model  4. 

The  vector  plot  for  model  3  shows  a  smooth  flow  of  air  to  the  engine  and  augmenter 
tube.  The  counterflows  first  observable  in  Figure  34  (  model  4  )  become  more  pro¬ 
nounced  in  Figure  35  (  model  5&.6  ).  It  appears  that  in  the  test  bay  a  higher  magnitude 
of  turbulent  viscosity  results  in  a  more  distortion  free  flow  path  to  the  engine  and 
augmenter  tube.  Velocity  vector  plots  at  the  augmenter  tube  exhaust  for  model  3  and 
model  4  are  shown  in  Figure  31  on  page  62  and  Figure  32  on  page  63.  The  velocity 
vector  plots  for  model  5&6  and  model  7  are  not  shown.  Model  5&6  is  identical  to  model 
4  both  in  magnitude  and  direction.  The  velocity  vectors  in  model  7  are  identical  in  di¬ 
rection  to  model  4  but  have  a  magnitude  1.5  times  greater.  The  gas  flow  in  model  3  exits 
to  the  atmosphere  at  a  velocity  of  30  and  is  uniformly  distributed  across  the  exhaust 
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opening.  The  gas  flow  in  model  4  exits  at  the  far  end  of  the  exhaust  opening  at  a  ve¬ 
locity  of  50  In  addition,  a  counterclockwise  vortex  is  present  over  the  exhaust 
opening.  As  in  the  test  bay,  a  higher  magnitude  of  turbulent  viscosity  results  in  a  more 
distortion  free  flow  field. 

The  temperature  contours  in  the  area  of  the  jet  engine  exhaust  for  model  3,  model 
4,  model  5&6,  and  model  7  are  shown  in  Figure  38  on  page  69  through  Figure  41  on 
page  72.  Temperature  contours  at  the  augmenter  tube  exhaust  for  model  3.  model  4. 
model  5&6,  and  model  7  arc  shown  in  Figure  42  on  page  73  through  Figure  45  on  page 
76.  The  temperature  contours  for  model  3  (vj-5:6.4)  shows  the  largest  temperature  de¬ 
crease  through  the  augmenter  tube  (  A  T  =  382  °  K  ).  The  temperature  decrease  through 
the  augmenter  tube  for  model  7  (Nr~4.4)  is  A  T  =  341  “  K.  for  model  4  )  is  A 

T  =  233  “  K  ,  and  for  model  5&6  (vr^0.06)  is  A  T  =  173  °  K.  The  latter  eight  figures 
show  that  the  temperature  decrease  through  the  augmenter  tube  increases  as  the  mag¬ 
nitude  of  turbulent  viscosity  increases. 

The  pressure  contours  in  the  area  of  the  jet  engine  exhaust  for  model  3.  model  4. 
model  5&.6  and  model  7  are  shown  in  Figure  46  on  page  77  through  Figure  49  on  page 
so.  The  pressure  contours  at  the  augmenter  tube  exhaust  for  model  3,  model  4.  model 
5&.6.  and  model  7  are  shown  in  Figure  50  on  page  81  through  Figure  53  on  page  84. 
The  pressure  contours  from  model  3  and  model  7  are  similar  in  that  the  pressure  fields 
change  from  negative  pressures  in  the  augmenter  tube  entrance  to  a  positive  pressure 
at  the  augmenter  tube  exhaust.  The  pressure  contours  for  model  4  and  model  5&6  are 
similar  in  that  the  pressure  fields  are  negative  from  the  augmenter  tube  entrance  to  the 
exhaust  stack.  It  appears  that  there  is  a  crossover  point  of  the  magnitude  of  the  turbu¬ 
lent  viscosity  in  regards  to  distribution  of  the  pressure  field.  For  large  magnitudes  of  \v 
(4.4  -  6.4  -^)  the  pressure  changes  from  negative  to  positive  values.  For  small  magni¬ 
tudes  of  Vj  (0.01  -  1.1  )  the  pressure  field  is  always  negative. 
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Figure  12.  Augmeiiter  Tube  Length  Study  Pressure  Profile 
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Figure  13.  Augincnter  Tube  Study  Temperature  Profile 
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Figure  15.  Augmeiiter  Tulie  Study  Pres 
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Figure  18.  Velocity  Vector  Plot  for  Model  1  (test  bay) 
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Figure  19.  Velocity  Vector  Plots  for  model  2  (test  bay) 
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Figure  22.  Tempemture  Contour  for  Model  1  (augineiiter  tube) 
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Figure  26.  Pressure  Contour  for  Model  1  (augnienter  lube) 
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Figure  30.  Turbulence  Study  Pressure  Profiles 
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Figure  31.  Turbulence  Study  Temperature  Profiles 
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Figure  33,  Velocity  Vector  Plot  for  Model  3  (test  bay) 
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Figure  34.  Velocity  Vector  Plot  for  Model  4  (test  bay) 
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Figure  35.  Velocity  Vector  Plot  for 


;i  5  &  6  (test  bay) 


Figure  36.  Velocity  Vector  Plots  for  Model  3  (exhaust) 
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Figure  41.  Tenipcinture  Contours  for  [Model  7  (engine  exhaust) 


Figure  44.  Temperature  Contours  for  Model  5  &  6  (exhaust) 


Figure  45.  Temperature  Contours  for  Model  7  (exhaust) 
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Figure  46.  Pressure  Contours  for  Model  3  (engine  exhaust) 
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Figure  47.  Pressure  Contours  for  Model  4  (engine  exhaust) 
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Figure  49.  Pressure  Contours  for  Modei  7  (engine  exiiaust) 
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Figure  52.  Pressure  Contours  for  Model  5&6  (exhaust) 
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Figure  53.  Pressure  Contour  for  Model  7  (exhaust) 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

1.  Augmenter  Tube  Length  Study 

The  results  of  the  augmenter  tube  length  parametric  study  shows  no  significant 
changes  of  pressure,  temperature,  or  velocity  fields  when  comparing  the  models  with  a 
10"  o  dilference  in  augmenter  tube  length.  System  efficiency  based  on  the  augmentation 
ratio  for  the  shorter  augmenter  tube  model  is  2%  higher  than  in  the  original  model. 
Unfortunately,  it  is  likely  that  the  sound  level  at  the  augmenter  tube  exhaust  is  inveisely 
proportional  to  the  augmenter  tube  length.  That  is,  decreasing  augmenter  lube  length 
will  result  in  increasing  noise  levels.  Existing  jet  engine  test  facilities  barely  meet,  and 
in  many  instances  exceed,  acceptable  noise  level  standards  (S5  dba  at  200  feet  from  the 
test  cell  ).  Until  a  new  cost  etTccUve  acoustic  absorption  material  with  a  higher  coclli- 
cient  of  absorption  is  developed  or  a  noise  reduction  system  independent  of  augmenter 
tube  length  is  developed  the  design  of  the  Hush  House  incorporating  a  shorter 
augmenter  tube  would  be  contcrproductive. 

2.  Turbulence  Study 

'I  he  selection  of  a  turbulence  model  as  initially  expected  has  a  dramatic  elfect 
on  the  final  How  field  solution. .'\s  the  magnitude  of  the  turbulent  viscosity  increases,  the 
decrease  of  temperature  and  velocity  through  the  augmenter  tube  increase.  Halfway 
through  the  augmenter  tube  (Z  distance  =  39.7  meters)  model  (v;-:t6.4)  has  about  a  50",. 
reduction  in  temperature  and  about  a  90"  reduction  in  velocity  compared  with  model 
0  (\72:0.0l )  which  has  about  a  20"o  reduction  in  temperature  and  about  a  35"  n  reduction 
in  velocity.  'I  he  effects  of  turbulent  viscosity  on  the  pressure  distribution  through  the 
augmenter  tube  are  not  as  dear  cut  as  the  effects  of  turbulent  viscosits  on  temperature 


and  velocity.  The  magnitude  of  tlie  turbulent  viscosity  has  a  very  significant  efl'ect  on 
the  pressure  distribution  but  a  simple  correlation  between  the  magnit  ude  of  the  turbulent 
viscosity  and  the  pressure  distribution  cannot  be  made. 

It  is  interesting  to  note  that  in  general,  the  results  of  model  7  (k-  e  turbulence 
model)  fall  in  a  range  between  the  results  of  model  3  and  model  4.  The  magnitudes  of 
the  turbulent  viscosity  ENUT  in  the  arbitrilly  selected  cell  defined  by  lX-1,  n'-4,  and 
1Z=  20  for  model  3.  model  7.  model  4,  and  model  5  are  as  follows; 

ENL  T  for  model  3  =  6.4^ 

tn- 

ENL’T  for  model  7  =  4.4-=^ 

ni- 

ENLT  for  model  4  =  1.1-=^ 

m- 

ENL'T  for  model  5  =  0.06  ^ 

in- 

The  approximate  value  of  ENUT  =  4.4  for  model  7  can  be  arrived  at  with  the  multi¬ 
plying  factor  MF  in  equation  (3.10)  set  to  the  order  of  80,000.  That  is  the  k  c  model 
results  could  be  approximated  by  the  temperature  dependent  turbulence  model  by 
proper  setting  of  MF  in  equation  (3.10). 

This  turbulence  study  was  not  intended  to  determine  a  turbulence  model  that 
most  closely  simulated  the  actual  flow  field  in  the  Hush  House.  The  study  was  intended 
to  simply  show  how  different  turbulent  viscosity  models  effected  the  flow  field.  1  he 
study  shows  that  to  obtain  a  good  numerical  solution  using  PHOENTCS,  it  is  necessarv’ 
to  obtain  more  experimental  data  in  order  that  the  best  turbulence  model  is  incorporated 
in  the  PHOENICS  simulation. 

B.  RECOMMENDATIONS 

The  final  objective  is  to  develop  a  valid  computer  model  of  the  Naval  Air  Station 
Jacksonville.  Florida  Hush  House.  Two  major  stumbling  blocks  must  be  overcome  in 
order  to  attain  the  final  objective.  The  first  and  foremost  is  to  obtain  more  pressure. 
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temperature  and  velocity  data  from  the  augmenter  tube  of  the  actual  Hush  House.  At 
present,  temperature  and  velocity  data  is  available  for  only  one  location  in  the 
augmenter  tube.  The  number  of  locations  for  collecting  the  data  must  be  sufficient  to 
show  accurate  pressure,  temperature,  and  velocities  profiles  through  the  augmenter  tube. 
Once  the  data  is  obtained  the  PHOENTCS  generated  pressure,  temperature,  and  velocit> 
profiles  should  be  matched  to  the  actual  profiles  by  manipulating  the  magnitude  of  the 
multiply  factor  .MF  in  the  temperature  dependent  turbulence  model.  The  temperature 
dependent  turbulence  model  is  recommended  over  other  turbulence  models  for  future 
use  in  model  simulations  of  the  Hush  House  due  to  the  relative  ease  of  being  able  to 
match  model  results  to  actual  Hush  House  data. 
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APPENDIX  A.  "Ql"  FILE  (TURBULENCE  STUD^') 


9V>V*9ViVV'r'>VVr 


Vf'sVVnV*Vf^V*VnrtrVfVf*>V'A’^V*'A’Vc’?V^V7VVr*nViVVfVr'>V‘;VVf**‘)V'>V9V*Vf 


*  NAVAL  AIR  STATION  JACKSONVILLE  FLORIDA  HUSH  HOUSE  TEST 

*  FACILITY.  THESIS  PROJECT  OF  LT  ERIC  A.  NICOLAUS  USCG 

*  WITH  THESIS  INSTRUCTOR  PROF  DAVID  SALINAS. 


•VfVf*VcVcVc^V^VVf 


■VfVfVfVf*^VVfVcVc*Vf*T>’oVVf^V»VVfVr*'}'f'>V‘ 


Vc’VoViVVfVfVf^VTVVrVf^V^Vr*' 


* 

* 

* 

’^V^V 


*Vf^\’*Vr*5V**:ViVVryrVfVfVpVpVfyrVpV<*V?VfVrVfVrVfiVV«VV?TV^V‘)VVf^V*yfiV*VeVf>VTV^V^VyfV»"'5V*'>V*:V**TViV*Vr9ViVVrVf 


*  LIST  OF  VARIABLES 

* 


*  WJET:  VELOCITY  OF  ENGINE  EXHAUST  (METERS/SEC)  * 

*  MJET:  MASS  FLOW  RATE  OF  GAS  (KG/SECOND/SQ  METER)  * 

*  TJET;  JET  EXHAUST  GAS  TEMPERATURE  (DEGREES  KELVIN) 

*  PJET:  EXHAUST  PRESSURE  OF  JET  ENGINE  (NEWTONS/SQ  METER) 

*  RHOJET:  DENSITY  OF  THE  GAS  (KG/CUBIC  METER)  * 

*  HJET:  ENTHALPY  (JOULES/KG  OR  NEWTON*METERS/KG)  * 

*  KEJET:  KINETIC  ENERGY  (JOULES  OR  NEWTON -METERS)  * 

*  GASCON:  GAS  CONSTANT  ( JOULES/(KG-DEG  KELVIN))  * 

*  TAMB:  AMBIENT  AIR  TEMPERATURE  (DEGREES  KELVIN)  * 

*  HAMB;  AMBIENT  ENTHALPY  (JOULES/KG)  * 

*  KEINJ:  TURBULENT  KINETIC  ENERGY  (JOULES)  OF  JET  * 

*  CSUBP:  SP  HT  @  CONST  PRE  ( JOULES/(KG-DEG  KELVIN)) 

*  EPINJ:  DISSIPATION  RATE  OF  TURBULENT  KE  (JOULES)  OF  JET 

*  KEINA:  TURBULENT  KINETIC  ENERGY  (JOULES)  AMBIENT 

*  EPINA:  DISSIPATION  RATE  OF  TURBULENT  KE  (JOULES)  AMBIENT  * 
PI:  PRESSURE  (NEWTONS /SQUARE  METER) 

*  HI:  ENTHALPY  (JOULES/ KG) 

*  Ul:  VELOCITY  IN  X-DIRECTION  (METERS/SEC) 

*  VI:  VELOCITY  IN  Y-DIRECTION  (METERS/SEC) 

*  Wl:  VELOCITY  IN  Z-DIRECTION  (METERS/SEC) 

*  TMPl:  TEMP  OF  THE  FIRST  PHASE  (DEGREES  KELVIN)  * 

*  RHOl:  DENSITY  OF  THE  FIRST  PHASE  (KG/ CUBIC  METER)  * 

*  PRESSO:  REFERENCE  ATMOSPHERIC  PRESSURE  (NEWTONS/SQ  METER)  * 

*  TMPIA:  GROUND  REF  TEMP  USED  IN  GRND2  (DEG  KELVIN) 

*  TMPIB:  1/CSUBP  USED  IN  GRND2  EQU  (KG-DEG  KELVIN/ JOULES)  * 

*  DRHIDP:  DELTA  RHO/DELTA  PRESSURE  W/RESPECT  TO  HI  * 

*  FACT:  CORRECTION  OR  "FUDGE"  FACTOR  * 

*  ENUL:  LAMINAR  KINEMATIC  VISCOSITY  (METERS  SQ/SEC)  * 

*  ENUT:  TURBULENT  KINEMATIC  VISCOSITY  (METERS  SQ/SEC)  * 


•V?Vpy?yr>VVIrypyfV’oV>VyryrVr^V^Vyryr>V’?VVf*'>Vyr>VyrVrVrVpVpyryryryryryt'yp'>V* 


GROUP  1.  RUN  TITLE  AND  OTHER  PRELIMINARIES 
TEXT(NAS  JACKSONVILLE  HUSH  CHAMBER) 

REAL( WJET , MJET , TJET , PJET , RHOJET , HJET , KEJET , GASCON , CSUBP , TAMB , HAMB ) 
REAL( KEINJ, EPINJ, KEINA, EPINA) 

GROUP  2.  TRANSIENCE;  TIME-STEP  SPECIFICATION  ( DEFAULT=STEADY ) 

GROUP  3.  X-DlRECTlOi'.  GRID  SPECIFICATION 
NX=10 

XFRAC(1)=-1  ;  XFRAC(2)=0.  31 
XFRAC(3)=4  ;  XFRAC(4)=0.  5475 
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XFRAC(5)=3  ;  XFRAC(6)=2.0 
XFRAC(7)=2  ;  XFRAC(8)=1.4 


GROUP  4.  Y-DIRECTION  GRID  SPECIFICATION 
NY=12 

YFRAC(l)=-3  ;  YFRAC( 2)=0. 562 
YFRAC(3)=2  ;  YFRAC(4)=0. 314 
YFRAC(5)=3  ;  YFRAC( 6)=0. 562 
YFRAC(7)=4  ;  YFRAC(8)=0. 75 


GROUP  5.  Z 
NZ=40 

ZFRAC(1)=-1 

ZFRAC(3)=1 

ZFRAC(5)=1 

ZFRAC(7)=1 

ZFRAC(9)=2 

ZFRAC(11)=1 

ZFRAC(13)=1 

ZFRAC(15)=1 

ZFRAC(17)=2 

ZFRAC(19)=1 

ZFRAC(21)=2 

ZFRAC(23)=3 

ZFRAC(25)=2 

ZFRAC(27)=7 

ZFRAC(29)=1 

ZFRAC(31)=4 

ZFRAC(33)=4 

ZFRAC(35)=4 

ZFRAC(37)=1 


-DIRECTION  GRID  SPECIFICATION 

;  ZFRAC(2)=1.0 
;  ZFRAC(4)=1.7 
;  ZFRAC(6)=2.  0 
;  ZFRAC(8)  =  1.7 
;  ZFRAC( 10)=3.  0 
;  ZFRAC( 12)=1. 5 
;  ZFRAC(14)=1.  0 
;  ZFRAC(16)  =  1.  7 
;  ZFRAC( 181=1.0 
;  ZFRAC(20)=1.  7 
;  ZFRAC(22)  =  1.  0 
;  ZFRAC(24)=0. 566666 
;  ZFRAC(26)=0. 35 
;  ZFRAC(28)=3. 00 
;  ZFRAC(30)=2.  0 
;  ZFRAC(32)=0. 75 
;  ZFRAC(34)=0. 575 
;  ZFRAC(36)=0. 5 
;  ZFRAC(38)=3.  0 


GROUP  6.  BODY-FITTED  COORDINATES  OR  GRID  DISTORTION 


GROUP  7.  VARIABLES  STORED,  SOLVED  &  NAMED 
S0LVE(P1,H1,U1,V1,W1) 

S0LUTN(P1,Y,Y,Y,N,N,N) 

STORE ( TMP 1 ,RH01 ,ENUT) 

GROUP  8.  TERMS  (IN  DIFFERENTIAL  EQUATIONS)  &  DEVICES 

GROUP  9.  PROPERTIES  OF  THE  MEDIUM 
TAMB=295. 

TJET=700 

PJET=0. 

MJET=197. 72 
PRESS0=101325. 

RH01A=0. 

GASCON=286. 

RH01B=1.  /GASCON 

RHO JET= ( ( PRE  S  S  0+P JET ) *RHO 1 B ) /T JET 

KJET=MJET/RHOJET 

TMP1A=TINY 

***  tMPIB  is  EQUAL  TO  1.  DIVIDED  BY  C-SUB-P;  USED  IN  TMP1=GRND5 
CSUBP=1004. 

H.4MB=TAMB*CSUBP 


TMP1B=1. /CSUBP 
HJET=CSUBP*TJET 

**  THE  NEXT  2  STATEMENTS  INVOKE  COMPRESSIBLE  FLOW  IN  GREXl 
TMP1=GRND2 
RH01=GRND5 
DRH1DP=GRND5 
ENUT=GRND 

KEINJ=0.  5*(0. 005*WJET)**2  ;  EPINJ=0. 09*KEINJ**1. 5/(0.  05*0.  5) 
KEINA=0. 5*(. 005*1. 0)**2  ;  EPINA=0. 09*KEINA**1. 5/( 0. 1) 

GROUP  10.  INTER-PHASE-TRANSFER  PROCESSES  AND  PROPERTIES 


GROUP  11.  INITIALIZATION  OF  VARIABLE  OR  POROSITY  FIELDS 
FIINIT(H1)=HAMB 

FIINIT(TMP1)=TAMB  ;  FIINIT(DEN1)=1. 0 
FIINIT(KE)=KEINA  ;  FIINIT(EP)=EPINA 
FIINIT(P1)=0.  0 
FIINIT(W1)=0.  0 
FIINIT(V1)=0.  0 
FIINIT(U1)=0.  0 

PATCH(WINIT,INIVAL,1,5,1,8,15,30,] .1) 

INIT(WINIT,W1,0.  0,0.  2*WJET) 

INIT(WINIT,H1,0.  0,CSUBP*0.  2*TJET) 

INIT(WINIT,KE,0. 0,KEINJ*0.  2) 

INIT(WINIT,EP,0. 0,EPINJ*0. 2) 

***  SEE  GROUP  13  FOR  POROSITY  FIELDS 


GROUP  12.  CONVECTION  AND  DIFFUSION  ADJUSTMENTS 


GROUP  13.  BOUNDARY  CONDITIONS  AND  SPECIAL  SOURCES 

**DEFINE  SOURCES  A.\D  SINKS  FOR  JET  ENGINE** 
COKPOR(0. O.CELL, 1,1,4,5,12, 14) 

PATCH(JETIN, HIGH, 1,1, 4, 5, 11,11,1,1) 

COVALC  JETIN , PI ,FIXFLU , -MJET) 

COVAL( JETI N , H 1 , ONLYMS , SAME ) 

PATCHC JETOUT,LOW ,1,1,4,5,15,15,1,1) 

COVALC  JETOUT , P 1 , FIXFLU , MJET) 

COVALC JETOUT, HI , ONLYMS ,HJET) 

COVALC  JETOUT ,W1 , ONLYMS ,W JET) 

COVALC  JETOUT , KE , ONLYMS , KE INJ) 

COVALC JETOUT , EP , ONLYMS , EP I NJ ) 

**EAST  BOUNDARY  TO  ATMOSPHERE** 

PATCHC OUTSIDE, EAST, 10, 10, 1,12, 1,40, 1,1) 

COVALC OUTSIDE, PI, 0.  1,0.  0) 

CO VAL( OUTS I DE , H 1 , F I XVAL , HAMB ) 

**TOP  OF  HUSH  HOUSE  AND  TUBE  INCLUDING  TUBE  LIP** 
CONPORC  0. 0 , NORTH ,1,9,12,12,1,19) 

CO.NPORC  0.  0 , NORTH ,  1 , 5  ,  -8 ,  -8 , 18 , 27 ) 

CONPORC  0.0, NORTH , 1 , 5 , - 2 , - 2 , 35 , 35 ) 

CONPORC  0. 0 , NORTH , 1 , 5 , -4 , -4 , 36 , 36 ) 

CONPORC  0. 0 , NORTH ,  1 , 5 , -5  ,  -5 , 37 , 37 ) 

CONPORC  0.0, NORTH , 1 , 5 , -6 . -6 , 38 , 38 ) 

CONPORC 0. 0, NORTH, 1,5, -7, -7,39,39) 
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**TOP  OF  MODEL  TO  ATMOSPHERE** 

PATCK( TOPOuT , NORTH ,1.10,12,12,20,40.1,1) 
COVAL(TOPOUT,P1,0. 1,0. 0) 

COVAL( TOPOLT , H 1 , ONLYMS , HAMB ) 

PATCHCTOPP .NORTH, 10,10,12,12,1,19,1,1) 
COVALCTOPP,P1,0. 1,0.  0) 

COVAL( TOPP , H 1 , ONLYMS , HAMB ) 

**FLOOR  OF  HUSH  HOUSE,  TUBE  AND  GROUND** 

CONPORC  0.0, SOUTH , 1 , 10 , 1 , > 1 , 1 , 40 ) 

**FORE  AND  AFT  BULKHEADS  IN  HUSH  HOUSE  AND  TUBE** 
CONPORC 0. O.LOW, 1,9, 1,12, -1,-1) 

CONPORC  0. 0 , LOW, 6 , 9 , 1 , 12 , -20 , -20) 

CONPORC 0.  O.LOW, 1,5, 9, 12, -20, -20) 

CONPORC  0. 0 ,LOW , 1 , 5 , 1 , 2 , -35 , -35) 

CO.NPORCO.  O.LOW,  1,5, 3, 4, -36, -36) 

CONPORC  0.0, LOW, 1,5, 5, 5, -37, -37) 

CONPORC 0. O.LOW, 1,5, 6, 6, -38, -38) 

CONPORC 0. O.LOW, 1,5, 7, 7, -39, -39) 

CONPORC 0. 0, HIGH, 1,5, 8, 8, -39, -39) 

**VERY  FRONT  MODEL  STRIP  TO  ATMOSPHERE** 

PATCHC  FRONT , LOW , 10 , 10 , 1 , 12 , 1 , 1 , 1 , 1) 

COVALC FRONT, P 1,0.  1,0.  0) 

COVALC  FRONT , HI , ONLYMS , HAMB ) 

**THE  VERY  END  OF  OUR  MODEL  TO  ATMOSPHERE** 

PATCHC ENDEND, HIGH, 1,10, 1,12, 40, 40, 1,1) 

COVALC ENDEND.P 1,0.  1,0.  0) 

COVALC  ENDEND , H 1 , ONLYMS , HAMB ) 

**OUTSIDE  EAST  WALL  IN  HUSH  HOUSE** 

CONPORC  0. 0 .EAST ,-9,-9,9,12,1,19) 

CONPORC 0. 0, EAST, -9, -9, 1,8, 1,1) 

CONPORC 0. 0, EAST, -9, -9, 1,8, 3, 3) 

CONPORC 0. 0, EAST, -9, -9, 1,8, 5, 8) 

CONPORC  0.0, EAST ,-9,-9,1,8,10,11) 

CO.NPORC  0.  0 , EAST,  -9,-9,1,8,13,14) 

CONPORC  0. 0 .EAST ,-9,-9,1,8,18,19) 


**OUTSIDE  OF  TUBE  PLUS  TUBE  LIP  IN  HUSH  HOUSE** 
CONPORC 0. 0, EAST, -5, -5, 1,8, 18, 39) 

GROUP  14.  DOWNSTREAM  PRESSURE  FOR  PARAB=.  TRUE. 

GROUP  15.  TERMINATION  OF  SWEEPS 
FSWEEP=1  ;  LSWEEP=90 

GROUP  16.  TERMINATION  OF  ITERATIONS 

LITERCP1)=40 

GROUP  17.  UNDER-RELAXATION  DEVICES 
REALCDELT.FACT) 

DELT=50.  0/CNZ*0.  2*WJET) 


91 


FACT=10. 

FACT=1.  00 

FACT=0. 75 
FACT=3. 5 
FACT=5. 

FACT=0.  50 
FACT=0. 1 

RELAX(P1,LINRLX,0.  3) 

RELAX(  U1 ,  FALSDT ,  FACpi^DELT) 

RELAX( V 1 , FALSDT , FACT*DELT) 

RELAX( W 1 , FALSDT , FACT*DELT ) 

RELAXCHl , FALSDT, FACT^DELT) 

RELAX( KE , FALSDT , FACT*DELT) 

RE  LAX(  E  P ,  FALSDT ,  FAC’P'^DE  LT  ) 

GROUP  18.  LIMITS  ON  VARIABLES  OR  INCREMENTS  TO  THEM 
VARMIN(Ul)=-300. 

VARMAX(U1)=1000. 

VARMIN(Vl)=-300. 

VARMAX(V1)=1000. 

VARMIN(W1)=-100. 

VARMAX(W1)=1500. 

VARMIN(TMP1)=250 

VA?vMAX(TMPl)=TJET 

VARMIN(H1)=1000. 

VARMAXC  HI )=CSUBP*TJET 
VARMIN(Pl)=-3.  5E+05 
VARMAXC Pl)=  3.5E+05 
VARMINCRH01)=.  05 
VARMAXC RHOl) =5. 

VARMINCKE)=1. E-10 
VARMAXC KE)=1.E 10 
VARMINCEP)=l.E-08 
VARMINCVIST)=l.E-20 
VARMAXCEP)=1.  3*CKEJE'F-*1.  5) 


GROUP  19.  DATA  COMMUNICATED  BY  SATELLITE  TO  GROUND 
USEGRD=T 

GROUP  20.  PRELIMINARY  PRINT-OUT 

GROUP  21.  PRINT-OUT  OF  VARIABLES 

GROUP  22.  SPOT-VALUE  PRINT-OUT 

**MONITOR  THE  OUTPUT  AT  DESIGNATED  POINTS** 
IXMON=l  ;  IYM0N=4  ;  IZMON=16 
TSTSWP=1  ;  LUPR3=6 
ABSIZ=1.0  ;  0RSIZ=1.0 

GROUP  23.  FIELD  PRINT-OUT  AND  PLOT  CONTROL 
NUMCLS=10*NXPRIN 
ITABL=3  ;  NPLT=1  ;  NZPRIN=1 
IFLTL=LSWEEP 
NPRINT=3 
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GROUP  24.  DUMPS  FOR  RESTARTS 
SAVE=T 

RESTRT(P1,U1,V1,W1) 

RESTRT(ALL) 

STOP 

// 


APPENDIX  B.  PROGRAM  FOR  COPYING  MVS  DATA  SETS 


A.  COPYTAPE 


//RATNER  JOB  (2529,9999) , 'XHUSH' ,CLASS=A 
//  EXEC  PGM=IEBGENER 

//SYSPRINT  DD  SYSOUT=A 

//SYS IN  DD  DUMMY 

//SYSUTl  DD  DSN=MSS.  S1541.  PHOENICS.  C0NV7,DISP=SHR 

//SYSUT2  DD  DSN=MSS.  S1541.  PHOENICS. DF09G,DISP=(SHR) 

// 
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on  on  no  on  on  on  on  on 


APPENDIX  C.  "GROUND"  FILE  (GROUP  9) 


Section  5  shows  the  code  required  to  evaluate  "ENL  "1"  (eqation 

Q^V'>VVf'5VyrVfV<’V«’:fr*VfTVV:^VVfVrVfV«'’>V’sVVr*VfVfVr^VVrVfVrVfVrVfiVVf>V‘3VVr^VV?Vf*^VVfVf^VVfVf7V‘jVV?VrVf>VfVf^VVfVrVfVfsV*VfVf5V 

c 

C---  GROUP  9.  Properties  of  the  medium  (or  media) 

C 

C  The  sections  in  this  group  are  arranged  sequentially  in  their 
C  order  of  calling  from  EARTH.  Thus,  as  can  be  seen  from  below, 

C  the  temperature  sections  (10  and  11)  precede  the  density 

C  sections  (1  and  3);  so,  density  formulae  can  refer  to 

C  temperature  stores  already  set. 

9  GO  TO  (91,92,93,94,95,96,97,98,99,900,901,902,903) , ISC 

Q*VfVrVfVfVfVrVcVfVfVf')VVfVrVfVfV<’*Vf*VfVrVfVfVfVfVfVf*VfVfVr>VVfVfVfVn>V^V*^VVr'jVVfVf“>^V“'*Vr*)V'>V^VVfVrVfyrVfVfVrVfVfVf^VVf* 

900  CONTINUE 

C  *  . SECTION  10 . 

C  For  TMPl.  LE.  GRND - phase-1  temperature  Index  AUX(TEMPl) 

RETURN 

901  CO.NTINUE 

*  . SECTION  11  . 

For  T.MP2.  LE.  GRND - phase-2  temperature  Index  AUX(TEMP2) 

RETURN 

902  CONTINUE 

*  .  SECTION  12  . 

For  ELI.  LE.  GRND - phase-1  length  scale  Index  AUX(LENl) 

RETURN 

903  CONTINUE 

*  . SECTION  13  . 

For  EL2.  LE. GRND - phase-2  length  scale  Index  AUX(LEN2) 

RETURN 

91  CONTI.NUE 

’■  . SECTION  1  . 

For  RHOl.LE.GRND -  density  for  phase  1  Index  AUX(DENl). 

965  CO.NTINUE 
RETURN 

92  CONTINUE 

*  . SECTION  2  . . 

For  DRHIDP.  LE.GRND---  D( LN( DEN) )/DP  for  phase  1  (DIDP). 

RETURN 

93  CONTINUE 

*  . SECTION  3  . 

For  RH02.  LE.  GRND - density  for  phase  2  Index  AUX(DEN2). 

RETURN 

94  CONTINUE 

*  . SECTION  4 . 

For  DRH2DP.  LE. GRND---  D( LN( DEN) )/DP  for  phase  2  (D2DP). 

RETURN 

95  CONTINUE 

*  .  SECTION  5  . 

For  ENUT. LE. GRND -  reference  turbulent  kinematic  viscosity. 

CALL  FN5(AUX(VIST) ,AUX(TEMP1) ,4. 886,5. 432E-2,-2.  426E-5,7.  931E-9, 
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no  no 


1-1. 104E-12) 

CALL  F.\25(AUX(VIST).l.  UE-6) 

CALL  FN27(AUX(VIST) ,AUX(DEN1)) 

CALL  FN25(AUX(VIST).1000. 0) 

RETURN 

96  CONTINUE 

C  * . SECTION  6  . 

C  For  ENUL. LE. GRND -  reference  laminar  kinematic  viscosity. 

RETURN 

97  CONTINUE 

*  .  SECTION  7  . 

For  PRNDTLC  ).LE.GRND - laminar  PRA.NDTL  nos.  ,  or  diffusivity. 

RETU’RN 

98  CONTINUE 


*  . SECTION  8  . 

For  PHINT(  ).LE.GRND -  interface  value  of  first  phase(FIIl). 

RETURN 
99  CONTINUE 


C  * . 

C  For  PHINT( 
RETURN 

Q  *  Vr  Vr  -sV  *  V?  Vr  Vc  Vf  Vr  Vr  Vr  Vr  V? 


.  SECTION  9  . 

).LE.GRND -  interface  value  of  second  phase(FII2) 
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APPENDIX  D.  PROGRAMS  REQUIRED  FOR  IMPLEMENTING 

GROUND 


A.  COPYGRD 

//COPYGRD  JOB  ( 2529, 9999 ),' COPYTMVS ' ,CLASS=A 

//  EXEC  PGM=IEBGENER 

//SYSPRINT  DD  SYSOUT=A 

//SYS IN  DD  DUMMY 

//SYSUT2  DD  UNIT=SYSDA, 

//  DSN=MSS.  S1541.  PHOENICS.  SRC(GROUND) ,DISP=SHR 

//SYSUTl  DD  * 

PROGRAM  MAIN 

C  THIS  IS  THE  MAIN  PROGRAM  OF  EARTH 

C  FILE  NAME  GROUND.  FTN . 16  July  1986 

C 

C  (C)  COPYRIGHT  1984,  LAST  REVISION  1986. 

C  CONCENTRATION  HEAT  AND  MOMENTUM  LTD.  ALL  RIGHTS  RESERVED. 

C  This  subroutine  and  the  remainder  of  the  PHOENICS  code  are 
C  proprietary  software  owned  by  Concentration  Heat  and  Momentum 
C  Limited,  40  High  Street,  Wimbledon,  London  SW19  5AU,  England. 

C 

C 

C  PROGRAM  MAIN 

C 

C  1  The  following  two  COMMON'S,  which  appear  identically  in  the 
C  satellite  MAIN  program,  allow  up  to  25  dependent  variables  to 

C  be  solved  for  (or  their  storage  spaces  to  be  occupied  by 

C  other  variables,  such  as  density).  If  a  larger  number  is 
C  required,  the  25 ' s  should  be  replaced,  in  the  next  8  lines, 

C  by  the  required  larger  number;  and  the  100  in  COMMON/FOl/ 

C  should  be  replaced  by  4  times  the  required  number.  Numbers 
C  less  than  25  are  not  permitted. 

C 

C0MM0N/LGE1/L1(25)/LGE2/L2(25)/LGE3/L3(25)/LGE4/L4(25) 

1/LDB1/L5(25)/IDA1/I1(25)/IDA2/I2(25)/IDA3/I3(25)/IDA4/I4(25) 

1/IDA5/I5(25)/IDA6/I6(25)/GI1/I7(25)/GI2/I8(25)/HDA1/IH1(25) 

1/GH1/IH2(25)/RDA1/R1(25)/RDA2/R2(25)/RDA3/R3(25)/RDA4/R4(25) 

1/RDA5/R5(25)/RDA6/R6(25)/RDA7/R7(25)/RDA8/R8(25)/RDA9/R9(25) 

1/RDA10/R10(25)/RDA11/R11(25) 

1/GR1/R12(25)/GR2/R13(25)/GR3/R14(25)/GR4/R15(25) 

1/IPIP1/IP1(25)/HPIP2/IHP2(25)/RPIP1/RVAL(25)/LPIP1/LVAL(25) 

1/ IFPL/ IPL0( 25 ) /RFPL1/0RPRIN( 25 ) /RFPL2/0RMAX( 25 ) 

1/RFPL3/0RMIN( 25 ) /RFPL4/CELAV( 25 ) 

LOGICAL  L1,L2,L3,L4,L5,DBGFIL,LVAL 
CHARACTER*4  IHl , IH2 , IHP2 ,NSDA 
C 

C0MM0N/F01/I9(100) 

COMMON/DISC/DBGFIL 
EXTERNAL  WAYOUT 
C 


ono  cjoo  non  nn  ononnnon  nnnn  nnnn  non 


2  Set  dimensions  of  data-for-GROUND  arrays  here.  WARNING:  the 
corresponding  arrays  in  the  MAIN  program  of  the  satellite 
(see  SATLIT)  must  have  the  same  dimensions. 

COMMON/LGRND/LG( 20 ) / IGRND/ IG( 20 ) /RGRND/RG( 100 ) /CGRND/CG( 10 ) 
LOGICAL  LG 

CHARACTER*4  CG 

3  Set  dimensions  of  data-for-GREXl  arrays  here.  WARNING:  the 
corresponding  arrays  in  the  MAIN  program  of  the  satellite 
(see  SATLIT)  must  have  the  same  dimensions. 
COMMON/LSG/LSGD(20)/ISG/ISGD(20)/RSG/RSGD(100)/CSG/CSGD(10) 
LOGICAL  LSGD 

CHARACTER*4  CSGD 

4  Set  dimension  of  patch-name  array  here.  WARNING:  the  array 
NAMPAT  in  the  MAIN  program  of  the  satellite  must  have  the 
SAME  DIMENSION. 

C0MM0N7NPAT/NAMPAT( 500) 

CHARACTER*8  NAMPAT 

Declare  local  CHARACTER  variables. 

CHARACTER  NDUM4*4  ,NDUM6*6 ,1;DUM15*15 

5  The  numbers  in  the  next  two  statements  (which  must  be  ident¬ 
ical)  indicate  how  much  computer  memory  is  to  be  set  aside 
for  storing  the  main  and  auxiliary  variables.  The  user  may 
alter  them  if  he  wishes,  to  accord  with  the  number  of 

grid  nodes  and  dependent  variables  he  is  concerned  with. 
COMMON  F( 500000) 

NFDIM=500000 


6 


Logical-unit  numbers  and  file  names,  not  to  be  changed. 


DBGFIL=.  FALSE. 

CALL  DSCEAR(14,LUPR3, ' 
CALL  DSCEAR(6,LUDUM, ' 
CALL  DSCEAR(-10,LUSDA, 
CALL  DSCEAR(-14,LUPR1, 
CALL  DSCEAR(21,LUDST, ' 


' ,15,NDUM15,-11,16) 
,4,NDUM4,9,33) 

’ ,4,NSDA,0,0) 

' ,15,NDUM15,0,0) 

' ,4,NDUM4,9,33) 


User  may  here  change  message  transmitted  to  logical  unit 
LUPR3 

CALL  WRIT40( 'GROUND  STATION  IS  GROUND.  FTN  11  JULY  86  ') 
CALL  MAINl ( NFDIM , LUPRl , LUPR3 , LUSDA , NSDA) 

CALL  WAYOUT(O) 

STOP 

END 

SUBROUTINE  GROSTA 

(C)  COPYRIGHT  1984,  LAST  REVISION  1986. 

CONCENTRATION  HEAT  AND  MOMENTUM  LTD.  ALL  RIGHTS  RESERVED. 
INCLUDE  (SATEAR) 

INCLUDE  (GRDLOC) 

INCLUDE  (GRDEAR) 

.  .  .  This  subroutine  directs  control  to  the  GROUNDS  selected  by 
the  satellite  settings  of  USEGRX,  NAMGRD  &  USEGRD. 

Subroutine  GREXl  contains  much  standard  material,  eg. 
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options  for  fluid  properties,  several  turbulence  models, 
wall  functions,  etc. 

IF(USEGRX)  CALL  GREXl 

.  ...  ESTER  is  for  electrolytic-smelter  modelling  of  the  Hall-cell 
and  Soderberg  types  used  in  the  reduction  of  aluminium. 

IF(NAMGRD.  EQ.  'ESTR')  CALL  ESTRGR 

.  ...  SCRS  contains  the  simple-chemical-reaction-model  of 

combustion,  the  theoretical  basis  of  which  is  found  in  the 
book  "Combustion  &  Mass  Transfer"  by  D  B  Spalding  (1979) 

This  ground  also  contains  geometrical  features  of  a 
simplified  can  combustor. 

IFCNAMGRD. EQ.  'SCRS')  CALL  SCRSGR 

.  . . .  A  more  advanced  model  of  a  combustor  is  given  in  COMBGR. 

IFCNAMGRD. EQ.  'COMB')  CALL  COMBGR 

.  .  .  .  WJETGR  shows  how  to  represent  non-isotropic  effects  in  the 
turbulence  of  a  wall  jet. 

IFCNAMGRD.  EQ.  'WJET')  CALL  WJETGR 

.  .  .  .  TRACGR  contains  software  for  tracking  fluid  interfaces  by 

means  of  a  set  of  imaginary  particles  which  follow  the  motion 

IFCNAMGRD.  EQ.  'TRAC')  CALL  TRACGR 

.  , . .  PARTGR  is  used  to  solve  for  the  motion  of  particles  slipping 
relative  to  the  host  fluid.  A  spectrum  of  particle  sizes  can 
can  be  represented.  Each  particle  is  characterized  by  a  size, 
an  interphase  friction  coefficient,  an  evaporation  rate  &  a 
temperature. 

IFCNAMGRD. EQ. 'PART' )  CALL  PARTGR 

.  . . .  RADIGR  provides  the  coding  sequences  required  to  activate 
the  so-called  six-flux  radiation  model. 

IFCNAMGRD. EQ. 'RADI')  CALL  RADIGR 

.  .  .  .  G.AUSGR  provides  the  Gauss -Seidel  solver  as  an  alternative 

to  the  whole-field  linear  equation  solver  provided  in  EARTH. 

IFCNAMGRD.  EQ.  'GAUS')  CALL  GAUSGR 

.  .  .  .  NOZLGR  provides  initial  conditions  &  special  print  out  for 
a  convergent-divergent  nozzle  case  for  which  body-fitted 
coordinates  are  used. 

IFCNAMGRD.  EQ.  'NOZL')  CALL  NOZLGR 
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.  .  .  .  AEROGR  provides  inlet  boundary  conditions  &  initial  conditions 
for  a  one-half  C  grid  for  an  aerofoil. 

IFCNAMGRD.EQ.  'AERO')  CALL  AEROGR 


POLRGR  specifies  uniform  flow  boundary  conditions  into 
a  polar  domain  of  360  degree  extent. 

IF(NAMGRD.  EQ.  'POLR')  CALL  POLRGR 

BTSTGR  contains  the  sequenses  used  in  conjunction  with 
the  BFC  test  battery. 

IFCNAMGRD.EQ. 'BIST')  CALL  BTSTGR 

TESTGR  contains  test  battery  sequences  used  in  conjunction 
with  the  test-battery  SATLIT  subroutine,  TESTST. 

IFCNAMGRD.EQ.  'TEST')  CALL  TESTGR 

SPECGR  is  a  generic  "special"  GROUND  the  name  of  which  can 
be  used  by  anyone  for  their  own  purposes. 

IFCNAMGRD.EQ. 'SPEC')  CALL  SPECGR 


.  . . .  The  model  ground  is  for  the  insertion  of  new  user  sequences. 
IFCUSEGRD)  CALL  GROUND 


.  .  .  .  The  data  echo  is  now  called  at  the  preliminary  print  stage. 


IFCIGR. NE. 20)  RETURN 
IFC- NOT.  ECHO)  GO  TO  20 

CALL  DATPRNCY,Y,Y,Y,  Y,Y,Y,Y,  Y,Y,Y,N,  Y,Y,Y,Y, 
1  Y,Y,Y,Y,  Y,Y,Y,Y) 


RETURN 


20  CALL  DATPRNCY,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N,N) 


RETURN 


END 


SUBROUTINE  GROUND 

CC)  COPYRIGHT  1984,  LAST  REVISION  1986. 

CONCENTRATION  HEAT  AND  MOMENTUM  LTD.  ALL  RIGHTS  RESERVED. 
INCLUDE  CSATEAR) 

INCLUDE  CGRDLOC) 

INCLUDE  CGRDEAR) 

INTEGER  HIGH , OLD , AUX , SOUTH , NORTH , EAST , WEST 
LOGICAL  STORE, SOLVE, PR I NT 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  USER  SECTION  STARTS: 


1  Set  dimensions  of  data-for-GROUND  arrays  here.  WARNING:  the 
corresponding  arrays  in  the  MAIN  program  of  the  satellite 
and  EARTH  must  have  the  same  dimensions. 

COMMON/LGR.ND/LGC20)/IGRND/IGC20)/RGRND/RGC100)/CGRND/CGC10) 
LOGICAL  LG 
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CHARACTER*4  CG 

2  User  dimensions  own  arrays  here,  for  example: 

DIMENSION  UUH( 10,10) ,UUC( 10,10) ,UUX( 10,10) ,UUZ( 10) 

PARAMETERC  MY=20 , MX=1 ) 

DIMENSION  GV1(MY,MX),GW1(MY,MX),GTMP1(MY,MX) ,GRH01(MY,MX) 

D IMENS I ON  GH 1 ( MY , MX ) , GENUL( MY , MX) , GMACH ( MY , MX ) , GHC ( MY , MX ) 

D I  MENS  I  ON  GST(  MY ,  MX  )  ,  GCF  ( MY ,  ILX  )  ,  GQ(  MY ,  MX ) ,  GTAU  ( MY ,  MX ) 

COMMON  CP,GAMA,RAIR,PRL,RFAC 

3  User  places  his  data  statements  here,  for  example: 

DATA  NXDIM,NYDIM/10,10/ 

4  Index  functions  for  GROUND-EARTH  variable  references. 
LOW(I)=NPHI+I 

HIGH(I)=2*NPKI+I 

0LD(I)=3*NPHI+I 

IN(I)=4*NPHI+I 

STOREC I )=M0D( I SLN( I ) , 2 ) .  EQ.  0 
SOLVE(I)=MOD(ISLN(I) ,3).EQ.  0 
PRINT(I)=M0D(IPRN(I) ,2).  EQ.  0 
S0UTH(I)=-(KF(I)-1) 

N0RTH(I)=-(KF(I)+1) 

EAST(I)=-(KF(I)+NY) 

WEST(I)=-(KF(I)-NY) 

5  Insert  own  coding  below  as  desired,  guided  by  GREXl  examples. 
Note  that  the  satellite-to-GREXl  special  data  in  the  labelled 
COMMON'S  /RSG/,  /ISG/,  /LSG/  and  /CSG/  can  be  included  and 
used  below  but  the  user  must  check  GREXl  for  any  conflicting 
uses.  The  same  comment  applies  to  the  EARTH-spare  working 
arrays  EASPl,  EASP2,. . . . EASPIO.  If  the  call  to  GREXl  has  been 
deactivated  then  they  can  all  be  used  without  reservation. 

IXL=IABS(IXL) 

IF(IGR. EQ. 13)  GO  TO  13 
IF(IGR. EQ. 19)  GO  TO  19 

GO  TO  (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21, 
122,23,24),IGR 


GROUP  1.  Run  title  and  other  preliminaries 


1  GO  TO  (1001, 1002), ISC 

1001  CONTINUE 
RETURN 

1002  CONTINUE 
RETURN 


GROUP  2.  Transience;  time-step  specification 

2  CONTINUE 
RETURN 

c 


lui 


C -  GROUP  3.  X-direction  grid  specification 

C 

3  CONTINUE 
RETURN 

c 

C -  GROUP  4.  Y-direction  grid  specification 

C 

4  CONTINUE 
RETURN 

Q7V7V***Vr*7V7VVfVr7V7V^7VVf7V7V7V7VVC7V7V7V7VVr7Sf7V*1^7V7V7Vyr7V7V7fc’'/r7V'jV7t'A'*7V*';fr7V7V7V*7V7VVr>V7V*7V7V7V7V7V7V7V^>V 


c 

C -  GROUP  5.  Z-direction  grid  specification 

C 

5  CONTINUE 
RETURN 

Q7V7V7V7V^7V7V7V7V7V7V7V7V7V*7V7V7V7V7V7V7V7V7V7V7V7V7VlV^*7V7V7V7V7V7VlV'A'7VlV7VVr7V7t>V>V7V'^7V7V7V7V>V7V7V7V7V7V7V7V7V7V7V7V 


c 

C -  GROUP  6.  Body-fitted  coordinates  or  grid  distortion 

C 


6  CONTINUE 
RETURN 

p7V7V7V7V7V7V7VV?7V7VVcVf7V7Vlfr7V7V*7V7V7V7V7V7V7V‘jV7VVr7V7V7V7V7V7Hr7V7V7V>V^7V7V7V7V'>V7V7V'3V'>V7V7V7V7V>V7V7V*7Vyr7V7VVf7V7V7V7V 

c 

C---  GROUP  7.  Variables  stored,  solved  &  named 

7  CONTINUE 
RETURN 


(^7VVf7V7V7VVfVfVf7V7V7Vyf7V7V7V’?VVf7VVfVfV?Vf*Vfy<’7V7V7V7V7V7V7V7V7VynV7V**yf**yf*yf*7V7V7V7V7Sf7Wr**7V7V*Vf7V7VVf7V7V7V 

c 

C---  GROUP  8.  Terms  (in  differential  equations)  &  devices 
C 


8  GO  TO  (81,82,83,84,85,86,87,88,89,810,811,812,813,814,815) 
1,ISC 

81  CONTINUE 


C  *  . SECTION  1 . 

C  For  UlAD. LE.GRND -  phase  1  additional  velocity  (VELAD). 

RETURN 
82  CONTINUE 


C  *  . SECTION  2 . 

C  For  U2AD. LE. GRND -  phase  2  additional  velocity  (VELAD). 

RETURN 


83  CONTINX^E 

C  * . SECTION  3 . 

C  For  VIAD. LE.GRND phase  1  additional  velocity  (VELAD). 

RETURN 

84  CONTINUE 

C  * . SECTION  4 . 

C  For  V2AD. LE. GRND phase  2  additional  velocity  (VELAD). 

RETURN 

85  CONTINUE 


C  * . SECTION  5  . 

C  For  WlAD. LE.GRND -  phase  1  additional  velocity  (VELAD). 

RETURN 
86  CONTINUE 


C  *  . SECTION  6 . 

C  For  W2AD. LE. GRND -  phase  2  additional  velocity  (VELAD). 


102 


o  o  no 


RETURN 

87  CONTINUE 

C  * . SECTION  7 - VOLUMETRIC  SOURCE  FOR  GALA 

RETURN 

88  CONTINUE 

C  *  .  SECTION  8  ---  CONVECTION  FLUXES 

RETURN 

89  CONTINUE 

C  *  .  SECTION  9  ---  DIFFUSION  COEFFICIENTS 

RETURN 

810  CONTINUE 

C  *  .  SECTION  10  ---  CONVECTION  NEIGHBOURS 

RETURN 

811  CONTINUE 

C  *  .  SECTION  11  DIFFUSION  NEIGHBOURS 

RETURN 

812  CONTINUE 

C  *  .  SECTION  12  ---  LINEARISED  SOURCES 

RETURN 

813  CONTINUE 

C  *  .  SECTION  13  ---  CORRECTION  COEFFICIENTS 

RETURN 

814  CONTINUE 

C  * .  SECTION  14  ---  USER’S  SOLVER 

RETURN 

815  CONTINUE 

C  *  .  SECTION  15  ---  CHANGE  SOLUTION 

RETURN 


C  *  Make  all  other  group-8  changes  in  group  19. 


c 

C---  GROUP  9.  Properties  of  the  medium  (or  media) 

C 

C  The  sections  in  this  group  are  arranged  sequentially  in  their 
C  order  of  calling  from  EARTH.  Thus,  as  can  be  seen  from  below, 
C  the  temperature  sections  (10  and  11)  precede  the  density 

C  sections  (1  and  3);  so,  density  formulae  can  refer  to 

C  temperature  stores  already  set. 

9  GO  TO  (91, 92, 93, 94, 95, 96, 97, 98, 99, 900, 901, 902, 903), ISC 

900  CONTINUE 


*  . SECTION  10 . 

For  TMPl. LE. GRND -  phase-1  temperature  Index  AUX(TEMPl) 

RETURN 

901  CONTINUE 

*  . SECTION  11 . 

For  TMP2. LE. GRND- -  phase-2  temperature  Index  AUX(TEMP2) 

RETURN 

902  CONTINUE 

C  * . SECTION  12 . 

C  For  ELI.  LE.  GRND - phase-1  length  scale  Index  AUX(LENl) 

RETURN 

903  CONTINUE 

C  * . SECTION  13  . 

C  For  EL2.  LE. GRND - phase-2  length  scale  Index  AUX(LEN2) 

RETURN 
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91  CONTINUE 

C  * .  SECTION  1  . 

C  For  RHOl.LE.GRND -  density  for  phase  1  Index  AUX(DENl). 

965  CONTINUE 
RETURN 

92  CONTINUE 

*  . SECTION  2 . 

For  DRHIDP. LE.  GRND---  D(LN(DEN))/DP  for  phase  1  (DIDP). 

RETURN 

93  CONTINUE 

*  . SECTION  3  . 

For  RH02.  LE.  GR-ND - density  for  phase  2  Index  AUX(DEN2). 

RETURN 

94  CONTINUE 

*  . SECTION  4 . 

For  DRH2DP.  LE.  GRND---  D(LN(DEN))/DP  for  phase  2  (D2DP). 

RETURN 

95  CONTINUE 

*  . SECTION  5  . 

For  ENUT. LE.  GRND - reference  turbulent  kinematic  '^iscosity. 

CALL  FN5( AUX( VIST) ,AUX(TEMP1) ,4. 886,5. 432E-2 , -2. 426E-5 ,7. 931E-9 , 
1-1. 104E-12) 

CALL  FN25(AUX(VIST),1.  OE-6) 

CALL  FN27(AUX(VIST),AUX(DEN1)) 

CALL  FN25(AUX(VIST) ,20000. 0) 

RETURN 

96  CONTINUE 

*  . SECTION  6 . 

For  ENUL.  LE. GRND —  reference  laminar  kinematic  viscosity. 

RETURN 

97  CO.NTINUE 

*  . SECTION  7  . 

For  PRNDTLC  ).LE.GRND---  laminar  PRANDTL  nos.  ,  or  diffusivity. 

RETURN 

98  CONTI.NUE 

*  . SECTION  8  . 

For  PHINT(  ).LE.GRND -  interface  value  of  first  phase(FIIl). 

RETURN 

99  CONTINUE 

*  . SECTION  9  . 

For  PHINT(  ).LE.GRND -  interface  value  of  second  phase(FII2) 

RETURN 


--  GROUP  10.  Inter-phase-transfer  processes  and  properties 

10  GO  TO  (101,102,103,104) , ISC 

101  CONTINUE 

*  . SECTION  1  . 

For  CFIPS.  LE.  GRND---  inter-phase  friction  coeff.  AUX(INTFRC). 

RETURN 

102  CONTINUE 

*  . SECTION  2  . 

For  CMDOT.  EQ.  GRND-  inter-phase  mass  transfer  Index  AUX(INTMDT) 

RETURN 
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103  CONTINUE 

*  . SECTION  3 . 

For  CINT(  ).EQ.  GRND - phasel-to-interface  transfer 

coefficients  (COIl) 

RETURN 

104  CONTINUE 

*  . SECTION  4 . 

For  CINT(  ).EQ.  GRND - phase2-to-interfdce  transfer 

coefficients  (C0I2) 

RETURN 

**'>V'>V*-,V-'riV5V******Vr**A****5V***********************')V**Vf*****Vr'>V***ilr** 


---  GROUP  11.  Initialization  of  variable  or  porosity  fields 

11  CONTINUE 
RETURN 

5V****Vf*iV*********iV**iV******')V*'>V*********5V****'>V>V****5V**VriV***iV*iV**** 

---  GROUP  12.  Convection  and  diffusion  adjustments 

12  CONTINUE 
RETURN 


---  GROUP  13.  Boundary  conditions  and  special  sources 


13  CONTINUE 

GO  TO  (130,131,132,133,134,135,136,137,138,139,1310, 

11311, 1312, 1313, 1314, 1315, 1316, 1317, 1318, 1319, 1320, 1321), ISC 
130  CONTINUE 


C .  SECTION 

131  CONTINUE 

C .  SECTION 

RETURN 

132  CO.NTINUE 

C . SECTION 

RETURN 

133  CONTINUE 

C .  SECTION 

RETURN 

134  CO.NTINUE 

C . SECTION 

RETURN 

135  CONTINUE 

C .  SECTION 

RETURN 

136  CONTINUE 

C .  SECTION 

RETURN 

137  CO.NTINUE 

C . SECTION 

RETURN 

138  CONTINUE 

C .  SECTION 

RETURN 


1 

2 


3 


4 


5 


6 


7 


8 


9 


coefficient  =  GRND 
coefficient  =  GRNDl 

coefficient  =  GRND2 

coefficient  =  GRND3 

coefficient  =  GRND4 

coefficient  =  GRND5 

coefficient  =  GRND6 

coefficient  =  GRND7 

coefficient  =  GRND8 
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139  CONTINUE 


SECTION  10 


coefficient  =  GRND9 


C 


C 


c 


c 


c 


c 


c 


c 


c 


c 


c 


c 


c 


1310 

RETURN 

CONTINUE 

1311 

RETURN 

CONTINUE 

1351 

CONTINUE 

1312 

RETURN 

CONTINUE 

1313 

RETURN 

CONTINUE 

1314 

RETURN 

CONTINUE 

1315 

RETURN 

CONTINUE 

1316 

RETl’RN 

CONTINUE 

1317 

RETURN 

CO.NTINUE 

1318 

RETURN 

CONTINUE 

1319 

RETURN 

CONTINUE 

1320 

RETURN 

CONTI.NUE 

1321 

RETURN 

CONTINUE 

RETURN 

SECTION  11 

SECTION  12 

SECTION  13 

SECTION  14 

SECTION  15 

SECTION  16 

SECTION  17 

SECTION  18 

SECTION  19 

SECTION  20 

SECTION  21 

SECTION  22 


coefficient  =  GRNDIO 

-  value  =  GRND 

- value  =  GRNDl 

.  value  =  GRND2 

-  value  =  GRND3 

-  value  =  GRND4 

.  value  =  GRND5 

-  value  =  GRND6 

-  value  =  GRND 7 

-  value  =  GRND8 

-  value  =  GRND9 

■ . value  =  GRNDIO 


-  GROUP  14.  Downstream  pressure  for  PARAB=. TRUE. 


14  CONTINUE 
RETURN 


---  GROUP  15.  Termination  of  sweeps 
15  CONTINUE 

C  *  Make  changes  for  this  group  only  in  group  19. 

RETURN 

QyfyfyfyfyfVfVfyfyr*VfVfyfyfyfVfyryfVfyfyfyf**Vfyry<yrVr*ycyeVryf*yryryfyf*VfVfyrViryryfyr*VrVfyryfVfVfVfyfyfVryp^VfVfVf 

c 
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C---  GROUP  16.  Termination  of  iterations 
C 

16  CONTINUE 

C  *  Make  changes  for  this  group  only  in  group  19. 

RETURN 

Q*>VVrVfVf*VrVfVfVf*Vr**iV*iVVfVrVc’3VVlrVfVr*VcVf*Vr^VVrVr*'*VfVf'>V!fr**>VVf**Vr’>VVr^V7V*Vr*iV^V‘5VVfTV'>VVrVr*Vc'>V'3V 

0 

C -  GROUP  17.  Under-relaxation  devices 

C 

17  CONTINUE 

C  *  Make  changes  for  this  group  only  in  group  19. 

RETURN 

QVc*^V*VcVf*VfVcVrV?VcVfVf*Vr*V<’*Vf^VVf*VfVf*^V^V*iV*TVVoV‘/r'jViViViVTVVc'jV'>V*>V'i^‘sSrTVyfVfVf*^V*'>V‘»ViVyr9V>V^>V^iV 

c 

C---  GROUP  18.  Limits  on  variables  or  increments  to  them 
C 

18  CONTINUE 

C  *  Make  changes  for  this  group  only  in  group  19. 

RETUR.N 

(^^V*VrVr***VfVr*Vf:^->VVr*Vp*****9V*Vf*^V*>VVf****VfV(rVfVc*^Vr'i!lr***iSf'>SrTV**'>Vyf*jV*>V*>iVyrVr*^*^Vfr*»VVe’'jV‘>V 


c 

C---  GROUP  19.  Special  calls  to  GROUND  from  EARTH 


C 

19  GO  TO  (191, 192, 193, 194, 195, 196, 197, 198), ISC 

191  CONTINUE 

C  *  . SECTION  1 - START  OF  TIME  STEP. 

RETURN 

192  CONTINUE 

C  * . SECTION  2 - START  OF  SWEEP. 

RETURN 

193  CO.NTI.NUE 

C  * . SECTION  3 - START  OF  IZ  SLAB. 

RETURN 

194  CONTINUE 

C  * . SECTION  4 - START  OF  ITERATION. 

RETURN 

195  CONTINUE 

C  * . SECTION  5 - FINISH  OF  ITERATION. 

RETURN 

196  CONTINUE 

C  * . SECTION  6 - FINISH  OF  IZ  SLAB. 

C  OPEN  (UNIT  =  10,  FILE  =  'DATA') 

C  CALL  GETYX(V1,GV1,MY,MX) 

C  CALL  GETYX(W1,GW1,MY,MX) 

C  CALL  GETYX( AUX( TEMPI ) ,GTMP1,MY,MX) 

C  CALL  GETYX(AUX(DEN1) ,GRH01,MY,MX) 

C  CALL  GETYX(H1,GH1,MY,MX) 

C  CALL  GETYX(AUX(VISL),GENUL,MY,MX) 

C  DO  1000  1=1, MY 

C  GSONIC  =  (GAMA*RAIR*GTMP1(I,1))**.  5 

C  GVEL  =  (GV1(I,1)**2  +  GW1(I,1)**2)**  5 

C  IF((IZ  .EQ.  10)  .AND.  (lY  .  EQ.  1)  .AND. 

C  KISWEEP  .EQ.  LSWEEP))THEN 
C  WRITEdO,*)  'GV1=',  GV1(1,1) 

C  WRITEdO,*)  'GTMPl  =' ,  GTMP1(1,1) 

C  WRITEdO,*)  'GWl  =',  GW1(1,1) 
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C  END  IF 

CIOOO  GMACH(I,1)  =  GVEL/GSONIC 
C  CALL  SETYX(C1,GMACH,MY,HX) 

C  CALL  GETONE(AUX(TEMPl),TINF,MY,l) 

C  CALL  GETONE(Cl,MINF,MY,l) 

C  CALL  GETONE(Wl,WINF,MY,l) 

C  CALL  GETONE(AUX(DENl),RHOI,MY,l) 

C  GTR  =  TINF*(1  +  (RFAC*(GAMA  -  1. ) )*(MINF**2)/2) 

C  DO  1010  1=1, MY 

C  GQ(I,1)=(GRH01(I,1)*GENUL(I,1)/.  7)*( (GH1( 1 , 1) -CP*TOAL)/2.  5E-4) 
C  GHC(I,1)  =  GQ(I,1)/(GTR  -  TWAL) 

C  GST(I,1)  =  GHC(I,1)/(RH0I*WINF*CP) 

C  GTAU(I,1)  =  GRH01(I,l)*GENUL(I,l)*GWl(I,l)/2.  5E-4 
Cl 010  GCF(I,1)  =  2*GTAU(I,1)/(RH0I*WINF**2) 

C  CALL  SETYX(C2,GHC,MY,MX) 

C  CALL  SETYX(C3,GST,MY,MX) 

C  CALL  SETYX(C4,GCF,MY,MX) 

C  CLOSE  (UNIT  =  10) 

RETURN 

197  CONTINUE 

C  * . SECTION  7 - FINISH  OF  SWEEP. 

RETURN 

198  CONTINUE 

C  * . SECTION  8 - FINISH  OF  TIME  STEP. 

RETURN 

c 

C —  GROUP  20.  Preliminary  print-out 
C 

20  CONTINUE 
RETURN 

Qyf>V***V.--,'fV,'*V.".VVf*******VfVf***Vf*****Vf**V<"*Vf******Vf  ************•!'?**■>■?**** 

C 

C -  GROUP  21.  Print-out  of  variables 

C 

21  CONTINUE 

C  *  Make  changes  for  this  group  only  in  group  19. 

RETURN 

QVc9V*:VVfVfVpVfVr**Vr**Vr‘jVV?V(rVf:V9ViV*Vf*^VVfVpVf*Vf’>V>V*'A-*V<Tl!r*^V'»V^V>V>V>V>V‘>>yr'»VVr7VVf^VV'rTV*>VVc*.V^V'>V*Vr 

c 

C -  GROUP  22.  Spot-value  print-out 

22  CO.NTINUE 

C  *  Make  changes  for  this  group  only  in  group  19. 

RETURN 

C*********Vf********************Vr******************************** 

C 

C -  GROUP  23.  Field  print-out  and  plot  control 

23  CONTINUE 
RETURN 

c 

C---  GROUP  24.  Dumps  for  restarts 
C 


24  CONTINUE 
RETURN 


END 

/* 

// 


B.  COMPGRD 

//COMPGRD  JOB  (2529,9999) , ’NIEBURG’ ,CLASS=B,REGION=2000K 
//*MAIN  LINES=(99) 

//CL  EXEC  FORTVCL, 

//  FARM.  LKED=( 'LET,NCAL,SIZE=(5000K,64K) .NOXREF' ) 

//FORT.  SYSIN  DD  DSN=MSS.  S1541.  PHOENICS.  SRC(GROUND) ,DISP=SHR 
//FORT. SYSLIB  DD  DSN=MSS. PHOENICS. COMMONS. FORT, DISP=SHR 
//LKED. SYSLIB  DD 
//  DD 

//  DD 

//LKED.  SYSLMOD  DD  DSN=MSS.  S1541.  PHOENICS.  LOAD( GROUND) ,DISP=SHR 

// 


C.  LINKEAR 

//LINKEAR  JOB  ( 2529 ,9999) , '*CHAM' ,CLASS=B,REGION=7000K 
//’■-MAIN  LINES=(99) 

//LINK  EXEC  PGM=IEWL,PARM='XREF, LET, LIST, SIZE=(2000K,256K) ,AMODE=24' 
//SYS  PR  I  NT  DD  SYSOUT=’'' 

//SYSUTl  DD  UNIT=SYSDA,SPACE=(CYL,(1,1)) 

//SYSLIB  DD  DS.N=MSS.  S1541.  PHOENICS.  LOAD ,DISP=SHR ,MSVGP=PHOENICS 

//  DD  DS.N=MSS.  PHOENICS.  EARTHL.  LIB, DISP=SHR,MSVGP=PHOENICS 

//  DD  DSN=MSS.  PHOENICS.  RESLIB. LIB ,DISP=SHR ,NSVGP=PHOENICS 

//  DD  DSN=MSS.  PHOENICS.  SATLIB.  LIB, DISP=SHR,MSVGP=PHOENICS 

//  DD  DS.N=SYS1.  VLSKMLIB.DISP*SHR,UNIT=3350 

//  DD  DSN=SYS1.  VF0RTLIB.DISP=SHR,UNIT=3350 

//SYSLMOD  DD  DS.N=MSS.  S1541.  PHOENICS.  LOAD,DISP=SHR ,MSVGP=PHOENICS 

//SYSLIN  DD  * 

INCLUDE  SYSLIB( GROUND) 

INCLUDE  SYSLIB(GREXl) 

INCLUDE  SYSLIB(MACDEP84) 

ENTRY  MAIN 
NAME  EARTHP(R) 

// 


109 


APPENDIX  E.  PROGRAMS  REQUIRED  FOR  COMPACTING  FILES 


A.  PHOENXOl  LOADLIB 

//RATNER  JOB  (2529 ,9999) PHOENIX  LOAD  LIB* ,CLASS==A 

//*  STEP  1  RENAME  CURRENT  PHOENIX  LOAD  LIBRARY  TO  TEMPORARY  LIBRARY 

//STEPl  EXEC  PGM=IDCAMS 

//SYSPRINT  DD  SYSOUT=A 

//SYS IN  DD  * 

ALTER  MSS.  S1541.  PHOENICS. LOAD  - 
NEVNAMECMSS.  S1541.  PHOENICS. LOADTEMP) 

/* 

//*  STEP  2  ALLOCATE  SPACE  FOR  NEW  PHOENIX  LOAD  LIBRARY 
//STEP2  EXEC  PGM=IEFBR14 

//DDl  DD  UNTT=SYSDA,DISP=(NEW,CATLG, DELETE), 

//  DCB=(RECFM=U,BLKSIZE=13030), 

//  SPACE=(13030,(400,15,7)), 

//  DSN=MSS.  S1541.  PHOENICS. LOAD 

/* 


B.  PHOEN\02  LOADLIB 

// 

//RATNER  JOB  ( 2529 , 9999 ),' PHOENIX  LOAD  LIB',CLASS=A 
//*  STEP  3  COPY  RENAMED  OLD  LOAD  LIBRARY  TO  NEW  LOAD  LIBRARY 
//STEP3  EXEC  PGM=IEBCOPY 
//SYSPRINT  DD  SYSOLT=A 

//IN  DD  DSN=MSS.  S1541.  PHOENICS.  LOADTEMP, DISP=( OLD, KEEP) 

//OUT  DD  DSN=MSS.  S1541.  PHOENICS.  LOAD, DISP=( OLD, KEEP) 

//SYS IN  DD  * 

COPY  I.NDD=IN,OUTDD=OUT 

/* 

// 


C.  PHOENX03  LOADLIB 


//RATNER  JOB  (2529,9999),'  PHOENIX  LOAD  LIB' ,CLASS=A 
//*  STEP  4  DELETE  RENAMED  OLD  PHOENIX  LOAD  LIBRARY 
//STEP4  EXEC  PGM=IDCAMS 
//SYSPRINT  DD  SYSOUT=A 
//ZYSIN  DD  * 

DELETE  MSS.  S 1541.  PHOENICS.  LOADTEMP  NONVSAM  SCRATCH 

/* 

// 


no 


D.  PHOENXOl  SOURCELB 


//RATNER  JOB  ( 2529 , 9999 ),' PHOENIX  SOURCE  LIB' ,CLASS=A 

//*  STEP  1  RENAME  CURRENT  PHOENIX  SOURCE  LIBRARY  TO  TEMPORARY  LIBRARY 

//STEPl  EXEC  PGM=IDCAMS 
//SYSPRINT  DD  SYSOUT=A 
//SYSIN  DD  * 

ALTER  MSS. S1541. PHOENICS.  SRC  - 
NEWNAMECMSS. S1541.  PHOENICS.  SRCTEMP) 


/* 

//*  STEP  2  ALLOCATE  SPACE  FOR  NEW  PHOENIX  SOURCE  LIBRARY 
//STEP2  EXEC  PGM=IEFBR14 

//DDl  DD  UNIT=SYSDA,DISP=(NEW,CATLG,DELETE) , 

//  DCB=(RECF.M=FB,LRECL=80,BLKSIZE=800), 

//  SPACE=(23440,(200,30,7)) , 

//  DSN=MSS. S1541. PHOENICS.  SRC 

/* 

// 


E.  PHOE.N.\02  SOI  RCELB 


//RATNER  JOB  ( 2529 , 9999 ),' PHOENIX  SOURCE  LIB' ,CLASS=A 

//■•■••  STEP  3  COPY  RENAMED  OLD  SOURCE  LIBRARY  TO  NEW  SOURCE  LIBRARY 

//STEP3  E.XEC  PGM=IEBCOPY 
//SYSPRI.NT  DD  SYSOUT=A 

//IN  DD  DS.N=MSS.  S1541.  PHOENICS.  SRCTEMP, DISP=( OLD, KEEP) 

//OUT  DD  DSN=MSS.  S1541.  PHOENICS.  SRC, DISP=( OLD, KEEP) 

//SYSIN  DD  * 

COPY  INDD=IN,OUTDD=OUT 
I'¬ 
ll 


F.  PHOE.NX03  SOL'RCELB 


//RATNER  JOB  (2529,9999),'  PHOENIX  SOURCE  LIB' ,CLASS=A 
//*  STEP  4  DELETE  RENAMED  OLD  PHOENIX  SOURCE  LIBRARY 
//STEP4  EXEC  PG.M=IDCAMS 
//SYSPRINT  DD  SYSOUT=A 
//SYSIN  DD  * 

DELETE  MSS. S1541. PHOENICS. SRCTEMP  NONVSAM  SCRATCH 

/* 

// 


in 


APPENDIX  F.  PROGRAM  'TORM09" 


"FORM09''  is  required  to  format  the  PHOEXICS  output  file  "DF09"  to  a  file  that 


can  he  read  bv  the  graphics  program  "PllO'l'ON". 

//FORM09  JOB  (2529,9999), 'CHAM' ,CLASS=B 
//*MAIN  LINES=(20) 

//  EXEC  FORTVCLG 
//FORT. SYSIN  DD  * 

C  PROGRAM  FORM09 
C 

C  Increase  PARAMETER  settings  if  more  variables  or  a  larger  grid 
C  is  required. 

C 

PARAMETERC  NPHI P=80 ) 

C  NB  NXYMAX  must  be  >=  N2  and  NX*NY 
PARAMETERC  NXYMAX=2  0  000) 

CHARACTERS  NAME(NPHIP)  ,MESS(10) 

CHARACTER*15  FILIN, FILOUT 

CHARACTER*!  ANS 

CHARACTER* 10  ACCESS 

LOGICAL  LBUF( 10) ,STORE(NPHIP) , DEBUG 

DIMENSION  RBUF( NXYMAX) , IBUF( 10) 

C 

100  F0RMAT(1X,7I10) 

101  F0RMAT(6(1PE13. 6)) 

102  F0RMAT(1X,79L1) 

103  F0RMAT(1X,19A4) 

C 

C  Reset  DEBUG=. TRUE,  to  get  error  messages.. 

DEBUG=. TRUE. 

C 

C  Set  Logical  Units.  .  . 

LUIN=9 

LUuUl=10 

LUBUG=6 

C 

C  Set  File  Names.  .  . 

C  Name  of  unformatted  restart  file 
C  WRITE(6,*)  '  RESTART  FILENAME?' 

C  READ(5,’(A)')  FILIN 

C  Name  of  formatted  restart  file 
C  WRITE(6,*)  '  FORMATTED  FILENAME?' 

C  READ(5,'(A)')  FILOUT 

C 

C  Direct  access  or  sequential? 

C  WRITE(6,*)  '  DIRECT  ACCESS  (PHIDA)  FILE?  Y/N' 

C  READ(5,'(A1)')  ANS 

C  IFCANS.EQ.  'Y'.OR.  ANS.EQ.  'Y')  THEN 

C  Version  1.4  default 
C  ACCESS=' DIRECT' 
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on  noooonoonnoooo 


LENREC=256 

ELSE 

Version  1. 3 

ACCESS=' SEQUENTIAL' 

LENREC=100 
END  IF 

WRITE(6,*)  '  DEFAULT  LENREC?  (100  VI. 3,  256  VI. 4)  -  Y/N' 
READ(5,'(A)')  ANS 

IF(ANS.NE.  'Y'.AND.  ANS.NE.  'Y')  THEN 
WRITE(6,*)  '  LENREC?  (WORDS  NOT  BYTES)' 

READ(5/(I)')  LENREC 
END  IF 

NBYTES  is  the  number  of  bytes  per  single  precision  word. 
NBYTES=4 

Record  length  in  bytes. 

LREC=NBYTES*LENREC 


OPEN  input  file: 

OPEN ( UN IT=9,F0RM=' UNFORMATTED' , STATUS= ' OLD ' , 

1  ACCESS=ACCESS , IOSTAT=IERR ,ERR=90) 

IF(IERR. NE. 0)  GO  TO  90 
IF( ACCESS. EQ. 'DIRECT' )  THEN 
READ(LUIN,REC=1)  (MESS( I) , 1=1 , 10) 

ELSE 

READ(LUIN)  (MESS(I) ,1=1,10) 

ENDIF 

IF( DEBUG)  WRITE(LUBUG,*)  MESS 
IF(ACCESS. EQ. 'DIRECT' )  THEN 
READ( LUIN , REC=2 )  ( LBUF( I ) , 1=1 , 4 ) 

ELSE 

READ(LUIN)  (LBUF(I),I=1,4) 

ENDIF 

IF(DEBUG)  WRITE(LUBUG,*)  (LBUF( I) , 1=1 ,4) 

IF( ACCESS. EQ. 'DIRECT' )  THEN 

READ( LUIN , REC=3 )  NX , NY , NZ , NPHI , ( I BUF( I ) , 1=1 , 6 ) , LENREC 
ELSE 

READ( LUIN )  NX , NY , NZ , NPHI , ( IBUF( I ) , 1=1 , 6 ) , LENREC 
ENDIF 

NXNY=NX*NY 

IF(NXNY.  GT.  NXYMAX)  GOTO  92 
IF(NPHI.  GT.  NPHIP)  GO  TO  93 
IF(DEBUG)  WRITE(LUBUG,*)  NX,NY,NZ 

OPEN  output  file:  FORMATTED,  SEQUENTIAL  and  of  FIXED  RECORD  LENGTH. 
0PEN(UNIT=10,F0RM=' FORMATTED’ ,STATUS=’ OLD' , 

1  ACCESS=' SEQUENTIAL' , I0STAT=IERR,ERR=9 1) 

IF(IERR. NE. 0)  GO  TO  91 
WRITE(LU0UT,103)  (MESS( I ) , 1=1 , 10) 

WRITE(LU0UT,102)  ( LBUF( I ) , 1=1 ,4) 

WRITE(LUOUT,100)  NX ,NY,NZ ,NPHI , ( IBUF( I ) , 1=1 ,6) .LENREC 
IF(ACCESS.EQ. 'DIRECT')  THEN 
READ(LUIN,REC=4)  RINNER 
ELSE 

READ(LUIN)  RINNER 
ENDIF 

WRITE(LUOUT,101)  RINNER 
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IF(ACCESS.EQ.  'DIRECT')  THEN 
READC  LUIN , REC=5 )  ( NAME( I ) , 1=1 , NPHI ) 

ELSE 

READC LUIN)  (NAME(I),I=1,NPHI) 

ENDIF 

C  Set  C36. . .  names  for  variables  51  to  NPHI  for  version  1.4,  to  avoid 
C  control  characters  etc]  (Unless  already  stored) 

IFCACCESS. EQ.  'DIRECT' .  AND.  NPHI.  GT.  50)  THEN 
READC LUIN, REC=10)  ( STOREC I ) , 1=1 ,NPHI ) 

DO  6  1=51, NPHI 
J=I-15 

IF(. NOT. STOREC I) )  THEN 
IFCI.LT.  100)  THEN 

WRITE(NAME(I),'("C"  ,12,"  ")')  J 
ELSE 

WRITE(NAME(I),'("C"  ,13)')  J 
ENDIF 
ENDIF 
6  CONTINUE 
ENDIF 

WRITECLUOUT, 103)  (NAME( I) , 1=1 , NPHI) 

IF(DEBUG)  WRITECLUBUG,*)  (NAME( I ) , 1=1 ,NPHI ) 

IFCACCESS. EQ.  'DIRECT' )  THEN 
READC  LUIN , REC=6 )  ( RBUF( I ) , 1=1 , NX) 

ELSE 

READC LUIN)  (RBUF(I),I=1,NX) 

ENDIF 

WRITECLUOUT, 101)  (RBUF( I ) , 1=1 ,NX) 

IFCACCESS. EQ. 'DIRECT' )  THEN 
READC  LUIN , REC=7 )  ( RBUF( I ) , 1=1 , NY) 

ELSE 

READC LUIN)  (RBUF(I) ,I=1,NY) 

ENDIF 

WR I TE ( LUOUT , 1 0 1 )  ( RBUF( I ) , I = 1 , NY) 

IFCACCESS. EQ. 'DIRECT' )  THEN 
READC  LUIN ,REC=8 )  ( RBUF( I ) , 1=1 ,NZ) 

ELSE 

READC LUIN)  (RBUF(I) ,I=1,NZ) 

ENDIF 

WRITECLUOUT, 101)  (RBUF( I ) , 1=1 ,NZ) 

IFCACCESS. EQ. 'DIRECT' )  THEN 
READC  LUIN , REC=9 )  ( RBUF( I ) , 1=1 ,NZ) 

ELSE 

READC LUIN)  (RBUF(I),I=1,NZ) 

ENDIF 

WRITECLUOUT, 101)  (RBUF(I) ,I=1,NZ) 

IFCACCESS. EQ. 'DIRECT' )  THEN 

READC LUIN, REC=10)  (STORE(I), 1=1, NPHI) 

ELSE 

READC LUIN)  (STOREC I), 1=1, NPHI) 

ENDIF 

WRITECLUOUT, 102)  (STOREC I) ,1=1, NPHI) 

IF(DEBUG)  WRITECLUBUG, 102)  (STOREC I ), 1=1 , NPHI ) 

NVAR=0 

DO  5  1=1, NPHI 
I F ( STORE ( I ) )  NVAR=NVAR+ 1 
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5  CONTINUE 
IREC=10 
DO  10  K=1,NZ 
DO  20  I=1,NVAR 
IF( ACCESS.  EQ.  'DIRECT' )  THEN 
IREC=IREC+1 
NVREC=NXNY/ LENREC 

IF(MOD(NXNY,LENREC).  NE.  0)  NVREC=NVREC+1 
DO  30  IVREC=1,NVREC 
JF=(IVREC-1)*LENREC  +  1 
JL=IVREC*LENREC 
JL=MINO(JL,NXNY) 

READ( LUIN , REC=IREC )  ( RBUF( J) , J=JF , JL) 

30  CONTINUE 
ELSE 

READ(LUIN)  (RBUF(J),J=1,NXNY) 

ENDIF 

WRITE(LUOUT,101)  (RBUF( J) ,J=1,NXNY) 

20  CONTINUE 
10  CONTINUE 

IF( DEBUG)  WRITE(LUBUG,1900) 

1900  F0RMAT(1X,'  NORMAL  STOP  IN  PROGRAM') 

GO  TO  999 

C  Error  trapping. . , 

90  IF(DEBUG)  WRITE(LUBUG, 1901) 

1901  FORMATdX,'  ERROR  IN  OPENING  INPUT  FILE’) 

GO  TO  999 

91  IF(DEBUG)  WRITE(LUBUG,1911) 

1911  FORMATC IX, 'ERROR  IN  OPENING  OUTPUT  FILE') 

GO  TO  999 

92  IFCNXNY. GT.NXYMAX)  WRITE(LUBUG, 1921) 

I F  C  N 2 . GT.  NX YM AX )  WR I TE ( LUBUG , 1 9  2  3 ) 

1921  FORMATC IX, 'SETTING  OF  NX*NY  TOO  LARGE;  RESET  IN  FORM09  AND 
1,/, IX, 'RE-CREATE  THE  LOAD  MODULE') 

1923  FORMATC IX, 'SETTING  OF  NZ  TOO  LARGE;  RESET  IN  F0RM09  AND  ' 
1,/, IX, 'RE-CREATE  THE  LOAD  MODULE') 

GO  TO  999 

93  WRITEC LUBUG, 193) 

193  FORMATC IX, 'SETTING  OF  NPHI  TOO  LARGE;  RESET  IN  FORM09  AND  ' 
1,/, IX, 'RE -CREATE  THE  LOAD  MODULE') 

999  CONTINUE 
C  CLOSE  files  .  .  . 

CLOSE C  LUIN , STATUS=' KEEP'  ) 

CLOSE CLUOUT,STATUS=' KEEP' ) 

STOP 

END 

//GO.  FT09F001  DD  DSN=MSS. S1541. PHOENICS. AUG90,DISP=C0LD) 

//GO.  FTlOFOOl  DD  DSN=MSS.  S1541.  PHOENICS.  DF09C ,DISP=COLD) 

/* 

// 
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